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EXECUTIVE  SUMMARY 


A  study  funded  by  DARPA,  was  conducted  to  understand  the  mechanism  of  transport  in 
ion  channels  and  their  possible  application  as  single  molecule  detection  sensors.  This 
report  details  the  motivation  for  such  a  study,  the  various  approaches  used  and  the  results 
obtained  so  far.  In  this  project,  a  3-D  Poisson  Nemst  Planck  approach  using  PROPHET 
to  simulate  ion  channels  like  porin  was  demonstrated  for  the  first  time.  Techniques  and 
tools  were  developed  to  obtain  experimental  data  for  ion  channels  -  specifically  porin. 
Continuum  equations  that  predict  selectivity,  sensitivity,  and  gating  based  on  mean  field 
theories  have  been  demonstrated  for  two  channel  types  of  different  structure.  Channel 
types  of  same  structure  have  not  yet  been  compared  in  detail.  Multi-scale  simulations  that 
predict  selectivity  and  sensitivity  have  been  started  and  early  work  has  been  completed 
and  published.  This  is  the  first  time  that  multiscale  simulations  have  been  done  for  ion 
channels.  Gating  has  not  yet  been  analyzed.  Experiments  showing  how  well  porin  and  its 
mutants  are  described  by  continuum  theories  and  by  multi-scale  simulation  are  finished 
on  two  mutants  and  well  underway  on  two  others,  of  known  structure.  The  sensitivity  of 
porin  to  mutation  and  ion  binding  has  been  established  showing  its  feasibility  as  an  ion 
sensor. 

Major  Findings: 

•  Currents  have  been  measured  in  a  single  protein  molecule  of  known  structure  in  a 
wide  range  of  solutions.  Previous  measurements  were  made  in  an  ensemble  of 
channels  of  unknown  orientation  and  number,  preventing  quantitative  analysis. 

•  Currents  have  been  compared  among  continuum  PNP  models  in  three  different 
and  independent  versions  and  experiments.  This  is  a  tremendous  increase  in 
understanding  because  no  one  before  had  used  a  theory  in  which  electrical 
potential  was  computed  from  Maxwell’s  equations. 

•  Currents  have  been  compared  among  four  models  of  selectivity  and  the  currents 
observed  in  three  channel  types  (Ca++,  Na+,  and  Cl  :  K  is  in  progress)  in  a  wide 
range  of  conditions.  This  is  an  enonnous  increase  in  capability  since  no  physical 
theory  has  matched  more  than  one  or  two  of  these  observations. 

•  Simulations  have  shown  the  existence  of  a  discrete  water  switch  in  carbon 
nanotubes  which  is  being  sought  experimentally  right  now.  This  is  a  qualitative 
change  in  understanding. 

•  Simulations  have  shown  how  to  put  ions  into  carbon  nanotubes  which  is  being 
sought  experimentally  right  now.  Simulations  have  shown  that  chemical 
functionalization  of  the  mouth  of  carbon  nanotubes  can  give  it  properties  like 
selectivity  in  an  ion  channel.  This  is  a  new  observation  which  has  lot  of  potential 
for  further  studies. 

•  Simulations  have  shown  how  to  embed  molecular  dynamics  data  into  continuum 
models.  Such  an  approach  is  very  helpful  to  bridge  the  time  scale  and  length  scale 
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gap  between  Molecular  Dynamics  simulations  and  continuum  simulations.  This  is 
a  qualitative  change  in  understanding. 

•  Standard  code  of  computational  electronics  has  been  used  to  predict  ion 
movement  through  channels.  Electronics  and  ionics  stand  united,  for  the  first 
time,  as  far  as  we  know. 

•  Self  consistent  simulations  of  ion  movement  (in  which  the  electric  field  is 
computed  directly  from  Poisson’s  equation)  have  been  made  to  predict  ion 
movement  through  gramicidin  channels.  Transport  Monte  Carlo  has  been 
implemented  for  channels. 
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INTRODUCTION 


Ion  channels  are  powerful  device  components  with  numerous  biological  functions.  They 
are  the  molecular  basis  for  the  movement  of  electrical  currents  across  membranes, 
resulting  in  signal  transduction.  Detecting  pathogens,  for  instance,  means  selecting  some 
unique  feature  associated  with  them,  and  then  creating  an  easily  measurable  signal,  such 
as  an  electrical  current,  of  that  detection.  The  key  features  of  ion  channels  that  give  rise 
to  their  extraordinary  range  of  biological  functions  include  variations  in:  1)  their 
selectivity  (what  ions  can  pass  through),  2)  their  conductance  (how  rapidly  can  ions  get 
through),  and  3)  their  sensitivity  (how  the  conductance  is  modulated  by  such  factors  as 
the  chemical  composition  of  their  environment,  the  transmembrane  voltage,  the 
membrane  surface  tension,  and  the  chemical  binding  of  the  ion,  etc.).  Some  of  the  various 
biological  functions  of  channels  that  potentially  have  enonnous  significance  for  device 
design  and  engineering  include: 

Signaling  and  computation— Ion  channels  are  the  current  carrying  and  regulatory 
molecules  of  the  nervous  system,  required  for  diverse  purposes,  including  sensory  input, 
cognition,  coordination  of  motion,  and  general  housekeeping  functions  of  the  body.  As 
elements  for  computation,  ion  channels  are  not  as  fast  as  solid  state  devices;  however, 
assemblies  of  ion  channels  as  organized  in  neuronal  membranes  are  capable  of  very 
complicated  and  specific  input-output  patterns.  This  presumably  is  the  basis  for  the 
brain's  excellence  in  pattern  recognition  and  in  the  coordination  of  complex  motions, 
capabilities  not  close  to  being  matched  by  human-constructed  computers  or  electronic 
devices. 

Triggers  for  cellular  events— Many  cellular  processes  are  triggered  by  electrical  signals 
generated  at  the  cell  surface  by  ion  channels.  This  includes,  but  is  not  limited  to,  muscle 
contraction,  glandular  secretion,  and  sperm  entry  into  the  egg. 

Electrical  power  generation— Some  species  of  fish  have  electric  organs  that  produce 
large  pulses  of  electrical  energy  which  stun  prey  in  the  nearby  water.  The  source  of  the 
electrical  current  and  voltage  is  ion  channels  in  membranes  that  are  stacked  in  series  with 
each  other,  much  as  cells  in  a  car  battery.  Selective  ion  channels  act  as  molecular 
batteries  that  transduce  the  chemical  potential  of  the  ionic  species  for  which  they  are 
selective  into  an  electromotive  force  that  can  drive  electrical  current. 

Energy  transduction— Proton  channels  are  critical  components  of  the  molecular 
apparatuses  for  transducing  energy  from  photons  (light)  into  a  form  that  can  be  used  for 
nutrient  synthesis  in  plants.  Proton  channels  are  equally  important  in  the  exchange  of 
energy  between  ATP  and  other  forms  of  useful  chemical  energy  in  the  cell.  The  ATP 
synthase  enzyme,  which  contains  in  one  protein  assembly  both  a  proton  channel  and 
catalyst  for  the  reaction  between  ATP,  ADP  and  free  P;,  is  almost  100%  efficient  and 
reversible.  Depending  on  the  proton  gradient  across  the  membrane  and  the  ATP/ADP/P; 
ratio,  it  can  use  the  proton  gradient  to  synthesize  ATP  or  use  an  ATP  surplus  to  pump 
protons  with  essentially  the  same  stoichiometry. 
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Fluid  pumping  and  filtration— The  universal  mechanism  for  moving  fluid  and 
electrolytes  across  boundaries  in  organisms  is  via  osmotic  pressure  created  by  regulated 
ion  movement  through  channels.  For  example,  channels  regulate  the  flow  of  water  across 
the  airway  epithelium  to  keep  it  moist,  and  across  the  intestinal  epithelium  to  regulate 
water  content  in  the  bowels.  Cystic  fibrosis  is  an  ion  channel  disease  which  results  in  a 
deficiency  in  the  movement  of  water  across  the  airway  epithelium,  and  cholera  is  a  lethal 
disease  caused  by  the  misregulation  of  ion  channels  in  the  intestinal  epithelium.  Another 
vital  example  of  this  aspect  of  channel  utility  is  the  function  of  the  kidney.  A  dramatic 
example  of  the  movement  of  fluid  via  osmotic  pressure  can  be  observed  in  root  cells 
which  provide  water  to  the  tops  of  the  highest  trees.  (Indeed,  it  seems  that  the  osmotic 
pressure  available  from  physiological  electrolyte  (about  300  milliosmolar)  is  a 
fundamental  limitation  of  tree  height.) 

Chemical  sensing— Ion  channel  proteins  often  contain  domains  that  are  selective  for 
various  chemical  species,  and  are  gated  by  binding  with  those  species.  For  example, 
chemically  sensitive  ion  channels  are  responsible  for  taste  sensing  in  the  tongue  and 
palate. 

Mechanotransduction— Ion  channels  are  responsible  for  mechanosensitivity  in  cells.  It 
has  been  shown  that  some  bacterial  channels  open  and  close  in  response  to  surface 
tension  changes  in  the  membrane  in  which  they  are  embedded.  In  other  cells,  the 
opening  and  closing  of  ion  channels  which  sense  mechanical  defonnation  are  mediated 
by  changes  in  cytoskeleton  tension,  rather  than  the  membrane  directly. 

In  summary,  ion  channels  in  nature  are  integral  components  of  a  myriad  of  biological 
processes  and  functions  that  have  clear  counterparts  in  human  technology.  Engineered 
ion  channels  are  thus  potentially  useful  components  in  micro-  and  nanotechnology 
applications. 

In  this  project  a  three-pronged  approach  to  build  the  foundations  of  ionic  channels  as 
engineered  devices  for  innovative  technological  applications  was  pursued.  A  team  of 
experts  to  address  the  computational,  experimental  and  engineering  (or  fabrication) 
aspects  of  ion  channels  was  assembled.  Specifically,  the  team  developed  some  of  the  first 
simulation  tools,  in  close  collaboration  with  the  experimental  effort,  for  multi-scale  and 
mixed-domain  (electrical,  chemical,  biological  and  microfluidic)  analysis  of  ion  channels 
to  obtain  fundamental  insights  into  the  selectivity,  conductivity,  and  sensitivity  of  ion 
channels  [19],  [6],  The  approach  is  summarized  in  Figure  1.  The  simulation  effort  relies 
heavily  on  the  experimental  and  fabrication  effort  to  extract  transport  coefficients  and  the 
structure  of  the  channel,  respectively.  Once  the  accuracy  of  the  transport  models  in  ion 
channels  has  been  established,  the  simulation  tools  will  enable  speculative  design  of 
miniaturized  devices  and  systems.  The  validation  of  the  simulation  enabled  speculative 
designs  will  be  established  by  the  experimental  and  engineering  effort. 
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Figure  1.  A  three-pronged  approach  to  enable  design  of  innovative  technology  based  on  ion 
channels 


As  shown  in  the  figure  above,  the  goal  of  our  project  is  to  demonstrate  Ion  Channels  as 
Biosensors.  This  goal  involves  several  tasks,  which  are  summarized  below: 

Task  1:  Continuum  simulators  for  ionic  channels  based  on  mean-field  theories 

1 . 1  Develop  fast  and  efficient  simulators  for  steady-state  analysis  of  continuum 
model  for  extraction  of  I-V  curves. 

1 .2  Create  reduced-order  models  for  fast  dynamic  analysis  of  the  continuum  model 

1 .3  Establish  accuracy  of  continuum  models 

1 .4  Calibrate  transport  coefficients  for  predictive  design 

Task  2:  Molecular  and  multiscale  modeling  of  ionic  channels 

2. 1  Develop  molecular  dynamics  (MD)  techniques  for  simulation  of  ionic  channels 

2.2  Conduct  Monte  Carlo  (MC)  and  combined  MC/MD  analysis  of  channels 

2.3  Extract  transport  coefficients  such  as  permittivity,  diffusion  coefficients  etc. 
from  molecular  analysis 

2.4  Conduct  multiscale  (combined  continuum/molecular)  modeling  of  ionic  channels 

Task  3:  Simulation  of  gating  phenomena  to  enable  simulation  driven  innovation  of 
chemical  and  biological  sensors 

3.1  Perform  continuum  modeling  of  various  gating  phenomena 

3.2  Conduct  molecular  modeling  of  various  gating  mechanisms. 
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3.3  Utilize  different  simulation  and  experimental  approaches  to  obtain  fundamental 
insights  into  selectivity  and  sensitivity  of  ion  channels. 

Task  4:  Experimental  effort  on  ionic  channels 

4.1  Obtain  I-V  data  for  gramicidin,  porin,  and  other  channels  to  enable  accuracy  of 
simulation  tools 

4.2  Obtain  experimental  data  on  transport  coefficients 

4.3  Perform  experiments  on  speculative  designs  enabled  by  simulation  tools 


The  following  pages  summarize  the  results  of  the  project.  The  first  part  of  the  report  is  on 
different  simulation  approaches  for  characterizing  ion  channels.  In  sections  3,  4  and  5 
results  from  Drift  Diffusion,  Monte  Carlo  and  Molecular  Dynamics  simulations  are 
discussed.  In  section  6  a  summary  of  the  theory  and  experimental  results  to  understand 
selectivity  is  given.  In  the  second  part  results  from  investigations  on  the  possibility  of 
incorporating  the  functionality  of  ion  channels  into  carbon  nanotubes  are  reported. 
Section  7  is  on  Molecular  dynamics  simulations  of  carbon  nanotubes  in  electrolytes  and 
Section  8  is  on  Quantum  correction  to  potentials  used  in  molecular  dynamics  simulations 
of  Carbon  nanotubes.  The  appendix  has  two  key  articles  published  due  to  this  project 
which  have  the  relevant  mathematical  equations. 
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CONTINUUM  THEORY:  DRIFT  DIFFUSION  SIMUFATIONS  OF 
PROTEIN  CHANNEFS 


Continuum  theories  are  based  on  a  flow  model  that  solves  for  the  transport  of  an 
equivalent  fluid  representing  the  macroscopic  properties  of  the  actual  particle  flow.  Self- 
consistent  continuum  models  include  continuity  equations  for  the  density  of  each  ionic 
species,  and  the  Poisson  equation  in  linear  or  nonlinear  form,  depending  on  the  approach 
used  for  the  solution.  The  flow  model  can  be  derived  starting  from  the  Langevin 
transport  equation  or  the  Boltzmann  transport  equation.  In  both  cases,  similar 
assumptions  must  be  made  of  carrier  distributions  near  equilibrium,  and  an  identical 
current  equation  with  a  drift  and  a  diffusive  term  is  obtained.  The  resulting  system  of 
equations  constitutes  the  Poisson-Nernst-Planck  model  of  transport,  which  is  the  ionic 
equivalent  of  the  drift-diffusion  Shockley  equations  used  in  semiconductor  devices  for 
electron  and  hole  transport.  The  continuum  approach  has  the  advantage  of  fast  solution 
methods,  and  good  flexibility  if  complicated  geometries  must  be  accommodated. 
Therefore,  it  is  very  suitable  to  explore  the  qualitative  behavior  of  structures,  by  varying  a 
vast  parameter  space,  fast  and  efficiently.  In  addition,  the  transport  model  is  based  on 
one  main  parameter  (mobility  or  diffusivity)  which  can  be  readily  calibrated  to  account 
for  higher  order  effects.  Continuum  models  are  also  necessary  to  formulate  compact  or 
reduced-order  models  for  the  study  of  arrays  of  elements.  For  structures  which  are  as 
small  as  ionic  channels,  however,  one  has  to  keep  in  mind  that  this  class  of  models 
provides  a  description  of  the  behavior  which  is  correct,  in  an  average  sense,  at  the 
tenninals.  The  carrier  densities  obtained  inside  the  structure  need  to  be  interpreted  as 
long-time  steady-state  averages  at  each  location,  and  are  not  representative  of  ionic 
charge  configurations,  which  might  be  established  at  any  one  instant  of  the  transport 
process.  The  same  consideration  holds  for  the  currents,  which  are  time-average  and  do 
not  resolve  the  fluctuations  exhibited  by  the  particle  flow.  Still,  continuum  models  are 
the  corner  stone  of  an  effective  hierarchical  simulation  strategy  and  can  be  efficiently 
calibrated  from  measurements  or  higher-order  models.  A  great  advantage  of  this  class  of 
approaches  is  that  complicated  phenomena,  like  ion  binding  and  release  on  the  wall  of  the 
pore,  can  be  represented  fairly  easily  in  terms  of  generation-recombination  terms  and  trap 
sites,  with  almost  no  consequence  in  the  overall  numerical  complexity. 

A  framework  for  three-dimensional  (3D)  drift-diffusion  simulations  of  ion  transport  in 
protein  channel  systems  has  been  developed,  based  on  the  simultaneous  solution  of 
Poisson’s  equation  [2],  [8].  The  detailed  derivation  of  the  PNP  equations  are  given  in  [17]. 
This  concept  describes  coulomb  interactions  between  all  mobile  and  fixed  charges,  and  a 
continuity  equation  for  each  ion  species,  describing  ion  permeation  down  an 
electrochemical  gradient.  Water,  protein  and  lipid  are  treated  as  continuum  regions 
characterized  by  specific  dielectric  constants.  Macroscopic  current  flow  is  resolved  by 
assigning  an  appropriate  mobility  and  diffusion  coefficient  to  each  ionic  species. 

The  drift-diffusion  model  is  implemented  using  PROPHET,  a  robust  commercial 
simulation  platform  based  on  state-of-the-art  numerical  engines.  The  correctness  of  the 
solvers  and  of  the  simulation  models  have  been  exhaustively  tested  and  validated  over  at 
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least  two  decades  by  industrial  and  university  groups.  The  non-linear  solution  procedure 
uses  a  robust  Newton  iteration,  and  discretization  of  the  flux  equations  uses  the  correct 
exponential  interpolation  realized  by  the  Scharfetter-Gummel  approximation. 

The  drift-diffusion  model  has  been  used  to  calculate  the  current-voltage  characteristics  of 
several  ion  channels,  including  OmpF  porin  and  some  of  its  mutants,  in  various 
electrolytes  [1]. 

The  PROPHET  simulator  has  been  used  to  compute  the  open  channel  current-voltage 
(IV)  curves  for  wild-type  porin  OmpF  and  its  mutant  G119D  in  a  range  of  solutions  of 
KC1  and  CaCF  [3],  [4],  It  is  straightforward  to  include  several  multi-valent  species  in  the 
PROPHET  framework  without  modifying  the  source-code.  For  the  purpose  of  matching 
current-voltage  measurements  at  the  system  level,  the  channel  is  treated  as  a  “black  box” 
and  a  spatially  uniform  diffusion  coefficient  is  used  as  the  fitting  parameter.  A  selection 
of  the  results  obtained  using  a  spatially  uniform  diffusion  coefficient  is  shown  in  figures 
2  and  3. 

The  mutation  of  wild-type  OmpF  to  G119D  increases  the  negative  charge  and  further 
narrows  the  cross-section  of  the  constriction  region  of  each  pore  of  the  trimer. 
Experimentally,  the  conductance  of  G119D  is  seen  to  be  15-40%  lower  than  the  wild- 
type.  PROPHET  was  used  to  study  the  effects  of  the  stearic  and  electrostatic  changes 
introduced  by  the  mutation,  separately.  Figure  4  shows  that  the  current-voltage  curves  of 
G119D  computed  both  with  and  without  the  extra  -2>\e\  charge  are  almost  identical.  This 
suggests  that  the  reduced  conductance  of  the  mutant  is  due  to  a  narrowing  of  the  pore 
constriction  rather  than  changes  in  the  distribution  of  pennanent  charge  that  accompany 
the  mutation. 


Figure  2.  Experimentally  measured  and  computed  current-voltage  curves  for  OmpF  in 
asymmetric  solutions  of  KC1:  0.25:0. 1M,  0.25:0.5M,  0.25:1M,  0.25:3M.  A  single  spatially 
uniform  diffusivity  D,  adjusted  to  give  the  best  fit  to  the  measured  curves,  was  used. 
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Figure  3.  Experimentally  measured  and  computed  current-voltage  curves  for  G119D  in 
50mMCaC12.  The  diffusivity  of  each  species  was  adjusted  to  give  the  best  fit  to  the  measured 
curves. 


Figure  4.  Comparison  of  drift-diffusion  model  with  measured  current-voltage  curves  for  G119D 
in  lOOmM  and  1M  KC1.  Measured  data  are  indicated  with  solid  (lOOmM)  and  dashed  (1M) 
lines,  the  simulation  data  are  indicated  with  symbols.  The  current-voltage  curve  obtained 
using  the  G119D  and  the  charge  distribution  of  wild-type  OmpF  (A)  is  almost  identical  to 
that  of  G119D  (■). 
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The  PROPHET  framework  for  specifying  an  arbitrary  space-dependent  diffusivity  for  all 
ion  species  is  now  complete.  Experimentally  known  ion  diffusivities  can  be  applied  in  the 
bulk-electrolyte  regions  and  coupled  with  values  that  vary  within  the  channel  pore, 
inferred  either  from  higher  order  physical  models  (e.g.,  Molecular  Dynamics)  or  from 
fitting  the  measured  current-voltage  curves.  Fitting  of  experimental  curves  using  a  more 
physical  model  of  diffusivity  is  inherently  more  difficult  due  to  the  restricted  tunability  of 
the  parameter  space. 


0  20  40  60  80 

Position  along  channel  (Angstroms) 

Figure  5.  Diffusion  coefficient  profiles  for  K+  (blue)  and  Cf  (red)  along  the  OmpF  porin  channel 
inferred  from  mean-square-deviation  calculations  of  Molecular  Dynamics  (symbols) 
together  with  a  polynomial  fit  (solid  lines)  and  3-level  average  fit  (dashed  lines). 
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Figure  6.  Current-voltage  curve  measured  for  OmpF  in  lOOmM  KC1  compared  with  current- 
voltage  curves  computed  using  (i)  two  diffusivity  profiles  shown  in  Figure  5  (polynomial  and 
3-level  average  fits  to  MD  data)  and  (ii)  uniform  diffusion  coefficient. 
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4  PARTICLE  METHODS:  MONTE-CARLO  SIMULATION  OF  ION 
CHANNELS 

The  accuracy  of  continuum  theories  is  questionable  for  nanoscale  systems  such  as  ion 
channels,  especially  when  experimental  data  is  not  available  and  diffusion  coefficients 
cannot  be  calibrated  easily.  To  address  this  issue,  we  have  developed  both  Monte  Carlo 
and  Molecular  Dynamics  approaches  for  ion  channels. 

Conduction  in  ion  channels  has  many  similarities  with  transport  in  ultra-scale 
semiconductor  devices  with  granular  doping.  We  have  adopted  Monte  Carlo  device 
simulation  strategies  to  pursue  the  study  of  ion  channels  from  a  device  engineering  point 
of  view.  This  was  accomplished  by  coupling  solution  reservoirs  to  the  conducting 
channel,  and  by  treating  the  movement  of  ions  self-consistently  in  3-D.  Our  previous 
work  indicated  that  main  problem  in  the  simulation  of  ionic  channel  conduction  is  that  for 
an  ion,  successful  traversal  of  the  channel  is  a  rare  event.  Many  attempts  of  ion  injection 
take  place,  but  the  vast  majority  of  particles  tend  to  return  back  to  the  original  solution 
reservoirs  at  the  ends  of  the  channel.  The  microscopic  particle  current  is  not  uniform  on 
short  time  intervals,  and  simulations  must  be  extended  to  millisecond  scale  to  approach 
the  averaging  that  takes  place  in  a  measurement  of  channel  conduction.  At  the  same 
time,  due  to  the  extremely  small  sizes  and  the  strong  Coulomb  forces,  the  time  step  must 
be  on  the  order  of  femtoseconds  for  a  correct  resolution  of  the  semi-classical  trajectories. 
Typical  semiconductor  device  models  require  tens  of  thousands  of  particles  but 
reasonable  results  are  obtained  in  several  picoseconds  of  simulation.  Therefore,  the  ionic 
channel  transport  problem  is  also  multiscale  in  time.  Monte-Carlo  computational  tools 
serve  an  important  function  by  bridging  the  gap  between  fast  continuum  simulations  and 
molecular  dynamics  simulations. 

A  Transport  Monte  Carlo  code  has  been  developed  to  simulate  ion  transport  in 
electrolytes  as  a  sequence  of  trajectories  interrupted  by  random  scattering  events.  Ion 
trajectories  evolve  in  three-dimensional  space  according  to  a  detailed  electric  field 
distribution,  obtained  self-consistently  by  solving  Poisson’s  equation  on  a  discrete  grid. 
Ion  positions  are  mapped  to  the  grid  to  generate  a  spatial  distribution  of  charge  density, 
which  provides  the  right-hand  side  of  the  discretized  Poisson’s  equation.  Discretization  of 
Poisson’s  equation  leads  to  a  truncation  of  short-range  electrostatic  forces.  Short-range 
charge-charge  interactions  are  included  explicitly,  by  evaluating  the  Coulomb  force 
within  a  sub-domain  surrounding  the  particles,  which  is  then  added  to  the  Poisson 
equation  force  (particle -mesh).  A  correction  term  to  eliminate  double  counting  in  the 
overlap  region  between  the  two  domains  is  also  required.  Ion  size  is  introduced  by 
including  a  Lennard-Jones  interaction  potential  for  ion-ion  interactions,  and  a  hard-sphere 
model  potential  for  impervious  boundaries.  This  creates  an  additional  force  that  prevents 
unphysical  ion  coalescence. 

Ion-water  electrostatic  interactions  are  handled  implicitly  by  assuming  the  water  to  be  a 
continuum  background  with  a  given  dielectric  permittivity.  Other  types  of  ion-water 
interactions  are  treated  stochastically  using  a  scattering  rate.  The  duration  between 
scattering  events  is  obtained  by  selecting  a  random  number,  uniformly  distributed  on  the 
unit  interval.  The  implicit  water  model  greatly  reduces  the  ensemble  size  compared  to 
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traditional  Molecular  Dynamics  simulations,  enabling  the  simulation  times  to  extend  to 
100s  of  nanoseconds. 


The  Transport  Monte  Carlo  code  is  currently  being  used  to  study  ion  transport  through 
the  gramicidin  channel.  Figure  7  shows  the  representation  of  the  gramicidin  channel  used 
in  the  3D  Transport  Monte  Carlo  model.  Electrolyte,  protein  and  lipid  are  represented  by 
the  blue,  green  and  red  regions  of  the  domain,  respectively.  Electrodes,  in  contact  with 
the  electrolyte,  are  placed  at  the  right  and  left  extremities  of  the  domain,  to  maintain  a 
fixed  bias  voltage  that  drives  the  ion  current.  The  trajectory  of  a  single  Na+  ion  trajectory 
crossing  the  channel,  under  a  250mV  bias,  is  shown.  Ion  trajectories  that  cross  the 
channel  are  extremely  rare  events.  In  this  5ns  simulation  in  1M  NaCl  only  two  Na+  ions 
were  observed  to  cross  the  channel.  No  anions  were  observed  to  enter  or  cross  the 
channel,  because  their  larger  ionic  radius,  which  is  comparable  to  the  average  channel 
radius.  When  ion  radius  is  ignored  both  Na  and  Cf  ions  can  traverse  the  channel, 
although  the  electrostatic  barriers  presented  by  the  particular  configuration  of  fixed 
charge  lining  the  pore  of  the  gramicidin  channel  limits  the  fraction  of  current  carried  by 
Cf. 


Figure  7.  Transport  Monte  Carlo  Simulation  of  NaCl  in  gramicidin 
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MOLECULAR  DYNAMICS  SIMULATIONS  OF  ION  CHANNELS 


The  finest-grained  level  of  description  feasible  for  describing  biomolecular  dynamics  is 
as  a  set  of  charged  balls  (the  atoms),  connected  by  springs  (the  bonds).  The  charges  on 
the  atoms,  the  polarizability  of  the  atoms,  the  sizes  of  the  atoms,  and  the  energetics  of 
deforming  the  bonds  are  collectively  known  as  the  force  field  of  the  system.  Force  fields 
for  biomolecular  simulations  are  derived  partly  from  spectroscopic  studies  on  small 
molecules  that  are  essentially  parts  of  macromolecules,  partly  from  quantum  mechanical 
calculations,  and  partly  are  adjusted  to  fit  known  constitutive  behavior,  such  as  self- 
diffusion  coefficient,  dielectric  constant,  etc.  Computing  molecular  dynamics  trajectories 
is  conceptually  simple,  since  atoms  simply  move  according  to  Newton’s  laws,  but  is 
computationally  quite  complex,  because  of  the  large  number  of  atoms  in  biomolecular 
systems  and  near  environment,  and  the  various  types  of  forces  involved.  There  are 
several  programs  that  are  conveniently  organized  to  simulate  biomolecular  systems. 
Three  of  the  most  commonly  used  are:  AMBER,  GROMOS,  and  CHARMM— our 
laboratory  uses  GROMOS. 

Given  a  structural  molecular  model  of  a  biomolecule,  the  essence  of  a  molecular 
dynamics  simulation  is  to  assemble  the  model  molecule  (i.e.,  assign  spatial  coordinates  to 
each  atom  in  the  molecule)  together  with  as  much  of  the  environment  as  is  needed  to 
provide  realistic  behavior,  initiate  thennal  motion  by  assigning  random  velocities  to  the 
atoms  with  a  Maxwellian  velocity  distribution  characteristic  of  the  temperature  at  which 
the  simulation  is  to  be  run,  and  then  let  the  system  move  according  to  Newton’s  laws.  By 
collecting  statistics  on  the  motions,  conformations,  and  the  energetics  of  the  system  as  it 
moves,  one  attempts  to  gain  insights  into  the  relationship  between  the  molecular  structure 
and  its  function. 

A  Molecular  Dynamics  study  of  ion  permeation  in  OmpF  porin  embeded  in  1  M  KC1 
Hydrated  Pahnitoyloleoylphosphatidylcholine  Bilayer  was  performed  [5], [6], [24].  The 
system  under  molecular  dynamics  (MD)  simulation  consists  of  a  trimer  of  OmpF,  633 
Pahnitoyloleoylphosphatidylcholine  (POPE)  lipid  molecules,  626  K+  cations,  592  Cl 
anions,  and  33283  water  molecules  (i.e.,  ca.  1  M  KC1).  The  ionization  states  (distribution 
of  charges  in  ionizable  residues)  in  porin  were  calculated  with  an  extended  theory  which 
includes  the  effect  of  dielectric  constant  on  the  shifts  of  pKa  values.  An  important  feature 
from  these  pKa  calculations  is  that  ASP312  is  deprotonated.  This  is  in  contrast  to  earlier 
model  (with  ASP3 12  protonated)  used  in  MD  simulation.  The  system  has  a  total  number 
of  144333  atoms.  MD  was  carried  out  by  the  particle-mesh-Ewald  method.  The 
simualtion  time  including  equilibration  is  6.5  ns.  Figure  8  shows  a  snapshot  of  the 
system. 

Diffusion  coefficients  as  a  function  of  the  z-position  (along  the  bilayer  nonnal)  for  K+, 
Cf,  and  water  were  calculated  from  the  MD  trajectories.  They  are,  respectively,  displayed 
in  Figures  9,  10,  and  11.  From  these  figures,  it  is  obvious  that  the  mobilities  of  the  ions  as 
well  as  water  inside  the  pore  are  smaller  than  those  outside  the  pore  region.  The  length  of 
the  pores  span  from  ca.  z  =  1.7  nm  (intracellular  side)  to  ca.  z  =  7.2  nm  (extra  cellular 
side).  Based  on  the  current  statistics,  the  z-profiles  of  diffusion  coefficients  and  ion 
densities  may  not  be  meaningful  to  draw  conclusions  on  ion  selectivity  and  penneation. 
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An  alternative  approach  is  to  apply  voltages  across  the  membrane  in  the  MD  simulations. 
A  preliminarily  investigation  of  the  feasibility  for  such  a  MD  simulation  was  completed. 
Movies  for  ion  translocations  across  the  channel  are  posted  at 
<http://peptide.ncsa.uiuc.edu/~schiu/porin.html> 


Figure  8.  A  snapshot  of  OmpF  in  POPE  from  MD  simulation 

In  the  future,  current-voltage  (I-V)  curves  will  be  calculated  from  MD  simulations  (with 
various  voltages  applied  across  the  membrane)  for  OmpF  porin  and  its  mutant  G119D. 
The  results  from  these  series  of  MD  simulations  may  give  a  better  insight  into  ion 
selectivity  and  permeation  of  porins. 
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Figure  9.  Diffusion  coefficient  profile  for  Cl"  ions. 


Figure  10.  Diffusion  coefficient  profile  for  K+  ion. 
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Figure  11.  Diffusion  coefficient  profiles  for  water. 
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THEORY,  SIMULATIONS  AND  MEASUREMENTS  OF  SELECTIVITY 


A  systematic  attempt  has  been  made  to  understand  selectivity  of  the  open  channel  using 
the  methods  of  the  physical  chemistry  of  concentrated  salt  solutions.  Computation  and 
analysis  of  a  system  as  complicated  as  an  ionic  channel  requires  tight  integration  between 
theory,  simulation,  and  experiment  [23].  If  theory  and  simulation  are  not  integrated 
closely  with  experiment  in  biophysics,  it  is  easy  to  do  irrelevant  work.  The  simulation 
and  experimental  approaches  were  validated  to  the  best  extent  possible. 

Theory  of  selectivity: 

The  first  contribution  was  to  use  the  simplest  model  of  concentrated  salt  solutions  that  fits 
data  (i.e.,  activity  coefficient)  over  the  entire  range  of  concentrations,  namely  the  Mean 
Spherical  Approximation  {MSA).  The  MSA  was  used  to  analyze  and  then  predict  behavior 
the  Z-type  Ca++  channel  [21].  This  channel  was  chosen  because  it  is  very  selective,  thus 
providing  a  severe  challenge  to  a  theory,  and  has  been  studied  in  considerable 
experimental  detail,  thus  providing  another  severe  challenge  to  a  theory.  The  channel  has 
also  been  well  characterized  by  molecular  genetics,  although  its  structure  is  not  known. 
Finally,  the  /.-type  Ca++  channel  is  among  the  most  important  channels  in  an  animal 
because  it  controls  the  heart  beat. 

The  MSA  theory  was  successful  in  fitting  a  large  fraction  of  the  experimental  data  on 
selectivity  of  the  Z-type  Ca++  channel  using  just  two  adjustable  parameters  (dielectric 
constant  and  channel  volume)  held  at  fixed  values. 

Monte  Carlo  (MC)  simulations  were  used  to  check  the  MSA  theory,  which  represent  the 
state  of  the  art  in  the  study  of  homogeneous  and  inhomogeneous  ionic  solutions,  to  check 
the  results  of  the  MSA.  The  results  checked.  Thus,  water  was  added  into  the  model,  using 
the  solvent  primitive  model  ( SPM)  to  compute  selectivity.  The  SPM  not  only  confirmed 
the  MSA  and  MC  results  but  also  was  able  to  fit  a  range  of  data  (e.g.,  Mg++)  that  had  not 
been  adequately  fit  previously  [16]. 

Work  was  begun  during  this  period  on  the  Density  Functional  Theory  of  selectivity,  as 
well  as  on  the  Na+  channel  and  the  CF  channel  [22], 

Measurements  of  selectivity: 

An  experimental  setup  has  been  completed  which  allows  measurement  of  single  channel 
currents  from  one  channel  molecule  in  a  wide  range  of  solutions.  The  apparatus  to  change 
solutions  on  both  sides  of  the  channel  has  allowed  13  solutions  changes  on  a  single 
molecule  in  the  best  situation.  It  is  nearly  100%  reliable  when  making  changes  on  just 
one  side  and  is  about  40%  reliable  when  changing  solutions  on  both  sides.  The  failure 
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mode  is  the  rupture  of  the  membrane,  which  is  believed  to  be  caused  by  electrical  and  not 
mechanical  transients. 

This  apparatus  allowed  speedy  measurements  of  the  current-voltage  (IV)  relations  of 
porin  in  some  1 5  KC1  solutions  of  various  concentrations  over  a  voltage  range  of  ±200 
mV.  This  data  set  is  now  the  mother  lode  of  experimental  data  of  IV  relations  on  channels 
that  must  be  fit  by  any  theory  of  penneation  or  selectivity  [21].  Experiments  were  started 
with  other  ions. 

Theory  of  Permeation: 

The  Poisson  Nemst  Planck  ( PNP )  was  extended  to  three  dimensions  using  the  spectral 
element  .  Several  papers  were  published  showing  that  a  direct  calculation  of  current 
voltage  (IV)  curves  was  possible  and  successful  [14],  [15],  [17],  [18],  [20].  This 

represents  the  first  use  of  spectral  elements  (and  perhaps  of  finite  elements)  to  a  channel 
problem. 

Nonequilibrium  Statistical  Mechanics: 

A  fundamental  problem  underlying  all  modeling  of  ion  channels,  and  much  modeling  of 
proteins,  is  the  correct  formulation  of  statistical  mechanics  of  nonequilibrium  systems. 
Work  has  been  started  on  a  formulation  of  statistical  mechanics  using  trajectories,  not 
states  as  the  fundamental  item  being  counted.  Subsequent  work  (after  the  agreement 
period)  has  shown  that  this  approach  promises  to  allow  a  general  application  of  the 
modem  theory  of  stochastic  processes  to  nonequilibrium  statistical  mechanics. 
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7 


ION  CHANNEL  BASED  CARBON  NANOTUBE  SENSORS 


Molecular  dynamics  (MD)  has  been  used  to  study  the  possibility  of  using  carbon 
nanotubes  as  “ion  channels”  i.e.  incorporating  the  functionality  of  ion  channels  in  carbon 
nanotubes  and  thereby  developing  nanopore  sensors  for  single  molecule  detection. 
Though  carbon  is  hydrophobic  MD  simulations  show  that  water  enters  the  nanotube.  The 
diffusion  rate  of  ions  into  the  tube  is  very  low  but  in  the  presence  of  an  external  electric 
field,  ions  enter  the  carbon  nanotube.  Additionally,  simulations  show  that  the  property  of 
selectivity  found  in  “ion  channels”  can  be  mimicked  by  chemical  functionalization  of  the 
ends  of  the  tube.  Currently  studies  are  being  done  to  simulate  these  under  biophysical 
conditions. 

Transport  of  molecules  through  macromolecular  pores  is  of  considerable  importance  in 
many  biological  and  nano-electromechanical  systems.  In  recent  years  engineered  ion 
channels  have  been  developed  to  function  as  single  molecule  detection  systems.  Figure 
12  shows  an  engineered  ion  channel  sensor.  In  an  applied  potential,  an  ionic  current  is 
carried  by  the  ions  that  bath  both  the  sides  of  the  lipid  bilayer.  When  the  target  molecule 
binds  to  the  binding  site  in  the  pore,  the  current  is  modulated  as  shown  in  the  trace.  The 
frequency  of  binding  reveals  the  concentration  of  the  analyte  and  the  duration  and 
amplitude  of  the  events  reveal  its  identity.  Though  they  have  significant  advantages 
including  high  sensitivity,  wide  dynamic  range  and  biocompatibility,  their  lack  of 
durability  makes  them  reliable  only  in  a  lab  setting. 


Figure  12.  (Left)  Schematic  of  single  molecule  detection  by  a  modified  alpha-Hemolysin 
channel:  Current  flow  decreases  on  analyte  binding.  (Right)  A  design  model  of  a 
functionalized  carbon  nanotube  in  a  slab  surrounded  by  water  molecules  and  K+  and  Cl- 
ions. 

A  review  of  recent  literature  showed  that  promising  options  for  such  a  device  that  could 
be  practically  realized  were  Gold  nanotubule  membranes,  ion  beam  etched  silicon  nitride 
membranes  or  single  carbon  nanotubes.  Carbon  nanotubes  being  rigid  and  having 
exciting  electrical  and  mechanical  properties  seem  ideal  for  an  integrated  circuit  chip 
design  but  fundamental  questions  have  to  be  addressed  regarding  transport  of  water  and 
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ions  through  it.  Quantitative  information  on  fundamental  modes  of  mass  transport  and 
transport  rates  have  been  obtained  through  carbon  nanotubes  of  diameters  of  about  150 
nm.  It  has  been  showed  by  molecular  dynamics  simulations  that  water  molecules  enter 
nanotubes  of  as  small  as  8.4  A  diameter  even  though  carbon  is  hydrophobic. 

This  project  has  recently  shown  for  the  first  time  the  possibility  of  transport  of  an 
electrolytic  solution  (KCL)  through  a  carbon  nanotube  by  molecular  dynamics 
simulations  [9].  A  short  13.4  A  long  uncapped  (16,  16)  single- walled  nanotube  was 
solvated  in  a  reservoir  of  1.85  M  KC1  and  simulated  using  GROMACS.  It  was  observed 
that  ion  occupancy  was  very  low  and  that  ions  do  not  diffuse  spontaneously  into  the  tube 
(Figure  13).  To  attract  the  ions  into  the  tube,  partial  charges  of  +/-0.38  e  were  placed  at 
atoms  on  the  rim  of  the  tube  to  create  a  dipole  -  the  positive  charges  being  on  top  rim  and 
the  negative  charges  on  the  bottom  rim.  The  tube  was  then  fixed  in  the  center  of  a  box  of 
length  33A  and  an  external  electric  field  was  applied  along  the  axial  direction  to  mimic 
the  membrane  potential.  It  is  observed  that  the  occupancy  increases  significantly  (See 
Figure  14).  Once  partial  charges  were  shown  to  increase  occupancy,  the  next  step  was  to 
replace  them  with  functional  groups. 


Gnripancy  nf  K-i-  inns  in  CNT 
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Figure  13.  Ion  occupancy  in  a  (16,16)  carbon  nanotube  (13.4  A  L,  21.696  A  dia)  fixed  in  a 
solution  of  1.85  M  KCL  Ionis  diffusion  is  very  low  without  electric  current  or  partial  charges. 

Functionalization  of  carbon  nanotubes  is  a  field  of  growing  interest.  Both  side  wall  and 
end  wall  functionalization  have  been  successfully  realized  in  recent  experiments.  In 
simulations,  an  asymmetric  functionalization  of  carboxylate  and  amino  residues  were 
used  instead  of  the  partial  charges  to  mimic  a  real  ion  channel  at  either  ends  respectively. 
The  partial  charges  were  obtained  from  GROMOS  force  field  database-  carboxylate 
group  from  glutamate  and  amino  group  from  lysine  residues  respectively.  The 
functionalized  carbon  nanotube  was  then  placed  in  a  slab  with  similar  properties  of  that 
of  a  lipid  bilayer  with  a  surrounding  bath  of  1.5  M  KCL  solution.  An  electric  field  of  0.15 
V/nm  was  used  to  drive  the  ions  through  the  tubes  and  the  trajectories  of  K+  and  Cl-  ions 
are  shown  in  Figure  15.  The  Chloride  occupancy  is  much  higher  than  potassium  ion 
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occupancy  even  though  the  net  partial  charges  at  the  mouths  are  same  in  magnitude.  The 
attractive  force  between  the  COO-  and  the  K+  ions  is  large  that  the  K+  ions  are  “bound” 
to  the  COO-  groups  at  the  mouth  reducing  the  K+  occupancy  in  the  tube.  Just  like  ionic 
channels,  the  functionalized  carbon  nanotube  also  seems  to  exhibit  selective  permeability 
of  anions  over  cations.  Further  investigation  is  needed  to  explain  the  phenomena  and  how 
it  can  be  utilized  to  control  the  rate  and  direction  of  ionic  current. 


Nn  nf  K+  ins  Hr  thfi  CNT 
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Figure  14.  Ion  occupancy  in  a  (16,16)  carbon  nanotube  (13.4  A  L,  21.696  A  dia)  fixed  in  a 
solution  of  1.85  M  KC1  with  external  electric  field  E=0.015  V/nm  and  partial  charges  of  +/- 
0.38  e  on  the  rim  atoms  (positive  on  the  top)  and  negative  on  the  bottom. 


K+  ton  liajectoiies  Cl- tan  trajectories 


Figure  15.  Motions  z(t)  of  individual  ions  (K+  on  the  left  and  Cl-  on  the  right)  shown  inside  the 
nanotube  along  the  nanotube  axis.  The  rate  of  chloride  ion  passage  is  higher  than  the  K+  ion 
passage  indicating  a  selectivity  of  anions  over  cations. 
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8  QUANTUM  CORRECTION  TO  VAN  DER  WAALS  INTERACTIONS 


Molecular  Dynamics  simulations  give  good  results  only  to  the  accuracy  with  which  the 
interatomic  interaction  potentials  are  correct[ll].  Since  at  molecular  scale,  quantum 
effects  also  might  be  important  along  with  classical  forces,  corrections  have  to  be  made 
to  the  various  interatomic  forces.  In  this  section,  we  present  a  quantum  correction  to  the 
van  der  Waals  interactions  which  can  be  used  to  model  nanoscale  devices. 

Many-body  corrections  to  van  der  Waals  interactions  from  semiclassical  Casimir  forces 
were  calculated  [10].  A  basic  analysis  of  the  role  of  van  der  Waals  terms  in 
electromechanical  systems  was  explored  and  its  significance  at  the  sub-nanometer  scale 
demonstrated  [13].  A  recent  theory  of  van  der  Waals  interaction  for  shells  of  pure  carbon 
is  based  on  universal  principles  fonnulated  in  1930’s.  A  new  approach  is  based  on  the 
quantum  electrodynamical  description  of  the  van  der  Waals/Casimir  forces[12].  A 
simple  and  effective  model  was  developed  to  estimate  the  many-body  contribution  due  to 
collective  modes  (plasmons).  This  contribution  is  believed  to  be  a  major  portion  of  the 
total  van  der  Waals  energy  because  of  the  high  oscillator  strength  of  the  plasmons.  The 
theory  reveals  many-body  terms  that  are  specific  for  various  low-dimensional  graphite 
nanostructures  and  are  not  taken  into  account  by  standard  one-body  calculations  within 
the  dispersionless  model  by  Lennard-Jones.  We  have  demonstrated  the  use  of  the  model 
for  several  systems  (shown  in  Fig.  16):  a  double -wall  nanotube  (A),  a  nanotube  on  the 
surface  (B)  and  a  pair  of  single  wall  tubes  (C).  A  significant  difference  has  been  shown 
for  the  dependence  of  the  van  der  Waals  energy  on  distance  (Fig.  17),  which  is  a 
consequence  of  quantum  correction  that  was  used. 

The  frequencies  of  the  plasmons,  which  are  mixed  by  the  Coulomb  interaction,  are 
explicitly  calculated  .  As  a  result  of  the  mixing,  the  total  system  energy  is  lowered  by  the 
van  der  Waals  contribution.  The  distance  dependence  of  the  new  many-body  correction 
tenn  has  a  fractional  exponent  5/2  for  the  tube-metal  cohesion  and  7/2  for  tube-tube 
interaction.  This  is  unlike  an  one-body  energy  given  by  LJ  (6-12)  potential.  It  was 


Figure  16.  Geometry  of  nanotube  systems  for  which  a  quantum  correction  to  van  der  Waals 
forces  has  been  calculated:  (A)  Double-wall  nanotube,  (B)  single  wall  nanotube  on  a  surface 
and  (C)  two  single  wall  nanotubes. 
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known  that  the  direct  summation  of  atom-atom  interactions  for  carbon  nanotubes  gives 
the  exponents  of  4  and  3  for  inter-tube  and  tube-substrate  cohesion,  respectively.  Thus, 
the  quantum  correction  term  is  a  totally  new  contribution  and  has  to  be  taken  into  account 
along  with  the  classical  term. 

Our  model  can  be  applied  to  simulations  of  biological  channels,  which  have  similar  one¬ 
dimensional  geometry.  This  allows  one  to  calculate  the  van  der  Waals  interactions  more 
accurately  which  in  turn  would  give  better  estimates  parameters  like  diffusion 
coefficients.  This  will  be  done  in  the  future. 


log  W,  eV 


Figure  17.  Van  der  Waals  potential  for  SWNT  on  a  metal  substrate  (upper)  and  between  two 
identical  tubes.  Nanotube  radius  was  0.7  nm  in  the  calculations. 
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CONCLUSIONS 


Numerical  and  experimental  investigations  have  been  performed  to  understand  and 
characterize  properties  of  ion  channels  in  order  to  use  them  as  biosensors.  This  project 
was  able  to: 

•  Demonstrate  use  of  mean  field  theories  to  predict  selectivity  and  sensitivity  based  on 
continuum  equations. 

•  Demonstrate  multi-scale  simulations  that  predict  selectivity,  sensitivity,  and  gating 
based  on  Langevin-Poisson  simulations. 

•  Perform  experiments  showing  how  well  porin  and  its  mutants  are  described  by 
continuum  theories  and  by  multi-scale  simulations. 

•  Show  biosensor  activity  of  porin  in  experiments. 

•  Demonstrate  carbon  nanotubes  can  be  used  to  mimic  ion  channels. 

•  Develop  quantum  correction  to  classical  theories  to  calculate  van  der  Waals  forces. 

The  results  obtained  through  this  report  have  led  to  breakthrough  advances  in  ion 
channels.  Specifically,  the  drawbacks  of  continuum  theory  have  been  pointed  out. 
Multiscale  methods  were  developed  to  overcome  these  drawbacks.  This  approach  is  the 
key  to  analyzing  the  conductance,  selectivity  and  sensitivity  functions  of  ion  channels. 
The  multiscale  methods  can  be  used  reliably  to  predict  conduction  through  a  variety  of 
ion  channels.  Finally,  the  development  of  nanopore  sensors  has  led  to  revolutionary 
ideas  for  detection  of  single  molecules.  The  results  obtained  can  be  used  in  developing 
engineered  ion  channels  or  artificial  ion  channel  like  nano-devices  in  single  molecule 
detection. 
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Hierarchical  Approach  to  Predicting  Permeation  in  Ion  Channels 
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ABSTRACT  A  hierarchical  computational  strategy  combining  molecular  modeling,  electrostatics  calculations,  molec¬ 
ular  dynamics,  and  Brownian  dynamics  simulations  is  developed  and  implemented  to  compute  electrophysiologically 
measurable  properties  of  the  KcsA  potassium  channel.  Models  for  a  series  of  channels  with  different  pore  sizes  are 
developed  from  the  known  x-ray  structure,  using  insights  into  the  gating  conformational  changes  as  suggested  by  a 
variety  of  published  experiments.  Information  on  the  pH  dependence  of  the  channel  gating  is  incorporated  into  the 
calculation  of  potential  profiles  for  K+  ions  inside  the  channel,  which  are  then  combined  with  K+  ion  mobilities  inside  the 
channel,  as  computed  by  molecular  dynamics  simulations,  to  provide  inputs  into  Brownian  dynamics  simulations  for 
computing  ion  fluxes.  The  open  model  structure  has  a  conductance  of  ~1 1 0  pS  under  symmetric  250  mM  K+  conditions, 
in  reasonable  agreement  with  experiments  for  the  largest  conducting  substate.  The  dimensions  of  this  channel  are 
consistent  with  electrophysiologically  determined  size  dependence  of  quaternary  ammonium  ion  blocking  from  the 
intracellular  end  of  this  channel  as  well  as  with  direct  structural  evidence  that  tetrabutylammonium  ions  can  enter  into 
the  interior  cavity  of  the  channel.  Realistic  values  of  Ussing  flux  ratio  exponents,  distribution  of  ions  within  the  channel, 
and  shapes  of  the  current-voltage  and  current-concentration  curves  are  obtained.  The  Brownian  dynamics  calculations 
suggest  passage  of  ions  through  the  selectivity  filter  proceeds  by  a  “knock-off”  mechanism  involving  three  ions,  as  has 
been  previously  inferred  from  functional  and  structural  studies  of  barium  ion  blocking.  These  results  suggest  that  the 
present  calculations  capture  the  essential  nature  of  K+  ion  permeation  in  the  KcsA  channel  and  provide  a  proof-of- 
concept  for  the  integrated  microscopic/mesoscopic  multitiered  approach  for  predicting  ion  channel  function  from 
structure,  which  can  be  applied  to  other  channel  structures. 


INTRODUCTION 

The  bacterial  potassium  channel  KcsA  marks  the  first 
permeation  pathway  high-resolution  structure  (Doyle  et 
al.,  1998)  from  the  superfamily  of  the  voltage-gated  ion 
channels,  which  may  be  present  in  most  living  cells  and 
underlie  innumerable  excitability,  transport,  signaling, 
and  osmoregulatory  functions  (Hille,  1992).  The  primary 
open-closed  gating  of  KcsA  is  by  pH  lowering  (Hegin- 
botham  et  al.,  1999)  with  a  voltage-dependent  distribu¬ 
tion  of  conducting  open  substates  (Meuser  et  ah,  1999). 
Perozo  et  ah  (1999)  used  electron  paramagnetic  reso¬ 
nance  spectroscopy  (EPR)  with  site-directed  spin  label¬ 
ing  to  show  that  during  gating  the  transmembrane  helix, 
which  forms  the  lumenal  surface  of  the  pore,  undergoes 
a  rotational  motion  that  provides  a  wider  opening  at  the 
intracellular  end  of  the  channel. 

A  number  of  computational  studies  using  methods  of 
molecular  dynamics  (MD)  (Allen  et  ah,  1999,  2000; 
Aqvist  and  Luzhkov,  2000;  Berneche  and  Roux,  2000; 
Shrivastava  et  ah,  2000;  Shrivastava  and  Sansom,  2000), 
Brownian  dynamics  (BD)  (Chung  et  ah,  1999;  Corry  et 
ah,  2000),  and  electrostatics  methodologies  (Roux  and 
MacKinnon,  1999;  Guidoni  et  ah,  2000)  have  been  done 
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on  KcsA,  providing  valuable  information  about  the  bio¬ 
physics  of  the  channel.  General  features  of  permeation  as 
shown  by  these  studies  include  the  stabilization  of  the 
selectivity  filter  by  ion  occupancy,  the  obligatory  single 
filing  of  ions  through  the  selectivity  filter,  and  the  pref¬ 
erential  attraction  of  cations  into  the  selectivity  filter  by 
a  combination  of  side  chain  and  backbone  charges  ex¬ 
posed  to  the  lumen.  It  has  also  been  calculated  that  ions 
are  significantly  attracted  into  the  wide  cavity  in  the 
center  of  the  channel  by  the  dipole  of  the  pore  helix 
(Roux  and  MacKinnon,  1999). 

However  there  has  not  yet  been  a  quantitatively  success¬ 
ful  description  of  ion  permeation  in  this  channel  that  pre¬ 
dicts  electrophysiological  properties  by  applying  physics 
with  no  arbitrary  parameters  directly  to  experimentally  mea¬ 
sured  structural  data.  The  purpose  of  the  present  study  is  to 
describe  the  implementation  of  a  hierarchical  strategy  (Ja- 
kobbson,  1998)  combining  multiple  computational  methods 
in  a  coordinated  way  in  which  all  the  parameters  in  the 
calculations  are  derived  directly  from  structural  data  with  no 
arbitrarily  adjustable  parameters  and  to  apply  this  strategy 
to  the  KcsA  potassium  channel.  As  a  high-resolution  open 
KcsA  channel  structure  is  not  yet  available,  we  have  con¬ 
structed  a  series  of  structures  of  different  opened  sizes  based 
on  our  best  interpretation  of  published  experimental  data  for 
how  the  open  channel  varies  from  the  x-ray  structural 
model,  which  is  functionally  closed.  The  resulting  com¬ 
puted  ion  fluxes  accurately  reproduce  known  electrophysi¬ 
ologically  measurable  properties  for  this  channel. 
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MATERIALS  AND  METHODS 
Molecular  modeling  of  the  open  channel 

Our  strategy  for  molecular  modeling  of  the  channel  takes  several  factors 
into  account.  1)  A  motion  of  the  Ml  transmembrane  helix  and  pore-lining 
M2  transmembrane  helix  tending  to  open  the  intracellular  end  of  the 
channel,  as  suggested  by  the  EPR  spectroscopy  of  Perozo  et  al.  (1999).  2) 
A  degree  of  opening  of  the  intracellular  end  of  the  channel  sufficient  to 
permit  tetrabutyl ammonium  ion  (TBA)  to  be  an  optimal  channel  blocker 
relative  to  other  quaternary  amines  (Meuser  et  al.,  1999;  Splitt  et  al.,  2000) 
and  to  enter  the  channel  completely  into  the  deep  cavity  (Zhou  et  al.,  2001). 
3)  The  preservation  of  the  shape  of  the  extracellular  vestibule,  as  suggested 
by  experimental  data,  that  differences  in  toxin  binding  affinity  of  open  and 
closed  K  channels  is  a  function  of  interactions  among  K  ions  in  the  channel 
rather  than  a  change  in  the  shape  of  the  vestibule  (Terlau  et  al.,  1999).  4) 
The  maintenance  of  hydrophobic  matching  of  the  transmembrane  helices 
with  the  suiTounding  implicit  membrane,  by  maintaining  a  constant  helix 
tilt  with  respect  to  the  bilayer  normal,  which  is  the  same  as  the  channel 
axis.  The  tendency  to  preserve  hydrophobic  matching  is  a  generally  ac¬ 
cepted  principle  of  protein-lipid  interactions,  albeit  not  always  perfectly 
obeyed  (Killian,  1998).  The  specific  modeling  process  is  described  below. 

The  x-ray  structure  of  KcsA  in  the  Protein  Data  Bank  is  an  L90C  mutant 
with  missing  residues  1  through  22  and  120  through  158  and  truncated  side 
chains  at  Arg-27,  Ile-60,  Arg-64,  Glu-71,  and  Arg-1 17  (Doyle  et  al.,  1998). 
This  structure  was  refined  by  completing  the  truncated  side  chains  and 
optimizing  their  positions  using  the  program  Modeler  (Sali  et  al.,  1990). 
The  ends  of  the  four  subunits  were  capped  with  the  nontitratable  groups 
CH3  — 0=0  and  NH — CH3  as  extensions  of  the  a-helices.  A  series  of 
energy  minimizations  was  done  using  a  distance-dependent  dielectric, 
leading  to  a  structure  that  satisfied  the  structural  criteria  of  Procheck 
(Laskowski  et  al.,  1993)  and  information-based  probability  density  func¬ 
tions  (Wall  et  al.,  1999). 

To  achieve  an  opening  of  the  intracellular  end  of  the  channel  consistent 
with  the  hydrophobic  matching  and  maintenance  of  the  external  vestibule 
shape,  the  extracellular  ends  of  the  transmembrane  helices  Ml  and  M2 
were  kept  fixed,  and  the  M2  helices  were  rotated  about  an  axis  parallel  to 
the  membrane  normal  (channel  axis).  The  Ml  helices  were  rotated  in  the 
same  fashion  as  necessary  to  eliminate  steric  overlap  with  M2.  The  Ml  was 
rotated  approximately  the  same  amount  as  was  the  M2. 

The  refined  crystal  structure  and  a  structure  formed  by  a  20°  rotation  of 
the  M2  were  inspected  for  gaps  in  packing  as  visualized  by  RasMol  (Sayle 
and  Milner- White,  1995)  in  space-filling  mode,  as  seen  in  Fig.  1,  A-D. 
When  viewed  from  the  extracellular  end  and  from  within  the  membrane 
plane,  both  the  closed  and  opened  structures  have  some  gaps  between  the 
Ml  and  M2,  which  are  presumably  filled  by  hydrocarbon  when  the  channel 
is  in  a  membrane.  Our  construction  procedure  may  widen  existing  gaps  but 
introduces  no  substantial  new  gaps  into  the  protein  structure. 

Fig.  1  E  shows  the  diameter  of  the  permeation  pathway  as  determined 
by  the  program  HOLE  (Smart  et  al.,  1993)  and  its  variation  along  the  axis 
in  the  Doyle  et  al.  (1998)  x-ray  structure  and  in  structures  of  various 
opened  sizes.  (For  convenience,  the  term  “opened  channel”  refers  to  the 
20°-opened  structure  unless  otherwise  stated.  The  20°  estimate  for  the 
degree  of  rotation  of  the  M2  helix  in  the  fully  opened  channel  is  based  on 
the  estimate  from  Perozo  et  al.,  1999.)  Two  very  narrow  regions  in  the 
x-ray  structure  are  found  at  the  selectivity  filter  and  near  the  intracellular 
end  of  the  channel.  The  intracellular  end  widens  to  an  average  diameter  of 
—  1 3.5  A  in  the  opened  channel.  This  size  is  consistent  with  experimental 
evidence  on  channel  blocking  by  quaternary  ammonium  ions  (Meuser  et 
al.,  1999;  Splitt  et  al.,  2000),  which  suggests  the  size  of  the  intracellular 
opening  for  the  maximally  conducting  substate  of  the  channel  is  larger  than 
TBA.  Further  evidence  that  the  channel  opens  widely  enough  to  permit 
TBA  to  enter  comes  from  a  recent  x-ray  structure  of  the  channel-TBA 
complex  showing  that  TBA  enters  completely  into  the  channel  (Zhou  et  al., 
2001).  (The  K  channel  structure  complexed  with  TBA  is  not  much  different 
from  the  closed  structure,  but  a  transient  opening  at  the  intracellular  end 


would  be  required  for  the  TBA  to  enter.)  Estimates  for  the  effective 
molecular  diameter  of  TBA  vary  from  a  minimum  of  ~9.5  A  (Splitt  et  al., 
2000)  to  a  value  of  11.6  A  based  on  an  idealized  geometry  (Huang  et  al., 
2000). 

It  should  be  noted  that  our  open  channel  structure  is  not  certain  in 
complete  detail.  However,  there  is  strong  experimental  evidence  that  the 
intracellular  end  of  the  channel  opens  to  approximately  the  extent  that  we 
suggest,  whereas  the  extracellular  vestibule  and  selectivity  filter  do  not 
change  in  dimension  or  conformation  during  channel  opening.  On  the  other 
hand,  the  M2  residues  are  generally  farther  away  from  the  channel  axis  in 
our  opened  channels  as  compared  with  the  closed  channel,  whereas  EPR 
experiments  (Perozo  et  al.,  1999)  suggest  that  some  of  the  residues  move 
closer  to  the  axis.  This  might  be  achieved  by  adding  to  the  motion  that  we 
have  suggested  a  rotation  of  the  helix  about  its  axis  and/or  a  kink  within  the 
helix  itself.  Given  the  absence  of  direct  evidence  for  specific  movements  of 
these  kinds  and  the  subtleties  inherent  in  interpretation  of  EPR  data,  we 
judged  it  not  useful  to  pursue  further  structural  modeling  without  more 
detailed  structural  data  on  the  open  state.  Because  of  the  residual  uncer¬ 
tainty  of  the  structure,  as  well  as  other  approximations  in  the  calculations 
that  we  discuss  later  in  the  paper,  we  emphasize  that  our  results  should  be 
considered  as  only  semiquantitative  representations  of  the  channel  behav¬ 
ior.  We  believe  the  major  conclusions  of  the  paper  depend  only  on  the 
intracellular  end  of  the  channel  opening  to  approximately  the  same  degree 
as  in  our  model  with  the  selectivity  filter  unchanged  and  will  therefore 
stand  as  long  as  the  correct  model  of  the  open  channel  has  these  features. 

Electrostatics:  pKa  calculations  and  potential 
profile  for  ion  permeation 

Electrostatics  calculations  were  used  to  calculate  the  ionization  states  of 
protein  residues  under  various  physiological  conditions.  We  used  the 
method  described  in  Gilson  (1993)  and  Antosiewicz  et  al.  (1994)  as 
implemented  in  the  University  of  Houston  Brownian  Dynamics  (UHBD) 
suite  of  programs  (Madura  et  al.,  1995)  used  to  calculate  electrostatic 
energies.  UHBD  requires  the  specification  of  dielectric  constants  for  the 
protein  and  surrounding  electrolyte.  To  emulate  the  low-dielectric  envi¬ 
ronment  around  the  channels  for  these  calculations,  a  layer  (thickness,  23.4 
A)  of  neutral  nonpolar  spheres  (radius,  2.0  A),  postulated  to  be  situated  in 
between  the  aromatic  residues  at  the  lipid- water  interfaces  (Cowan  et  al., 
1992),  was  placed  around  each  structure,  as  shown  in  Fig.  1  F.  The  spheres 
were  packed  at  hydrocarbon  densities  out  to  ~37  A  from  the  channel  axis. 
The  protein  dielectric  constant  was  set  to  20,  which  appears  to  be  the 
optimum  for  computing  the  ionization  states  of  residues  in  many  proteins 
(Antosiewicz  et  al.,  1994),  and  the  nonpolar  spheres  were  assigned  a  value 
of  20  as  well.  Outside  the  regions  delineated  by  protein  and  membrane  the 
dielectric  constant  was  set  to  80  to  represent  water.  UHBD  automatically 
accounts  for  the  image  charges  arising  from  introduction  of  a  dielectric 
interface. 

The  partial  charge  parameters  used  in  the  calculation  of  electrostatic 
energies  were  taken  from  the  CHARMm  (Brooks  et  al.,  1983)  parameter 
set  and  radii  from  the  OPLS  (optimized  parameters  for  liquid  systems) 
parameter  set  (Jorgensen  and  Tirado-Rives,  1988),  following  Antosiewicz 
et  al.  (1994).  Reference  pKa  values  corresponding  to  the  pKa  values  of 
model  compounds  in  solution  for  each  titrating  residue  type  are  taken  from 
Antosiewicz  et  al.  (1994).  A  modified  approach  (Tanford  and  Roxby, 
1972)  was  used  to  determine  effective  pKa  values  in  which  an  ionization 
polynomial  is  evaluated  exactly  for  clusters  of  residues  with  significant 
charge-charge  correlations  and  a  mean  field  approximation  is  used  for  less 
significant  intercluster  interactions  (Gibas  et  al.,  1997).  The  cluster  method 
requires  that  each  histidine  be  assigned  a  tautomer,  and  in  all  cases  the 
tautomer  in  which  Ne2  is  the  titrating  site  was  used.  The  calculations  used 
a  finite-difference  method  for  computing  electrostatic  potentials  at  the 
lattice  points  of  a  grid  of  spacing  1.2  A  surrounding  the  titratable  site. 
Focusing  grids  of  sizes  1.0,  0.75,  and  0.25  A  were  used  to  refine  the 
potentials.  (The  position  and  orientation  of  the  lattices  were  fixed  with 
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FIGURE  1  Packing  of  transmembrane  helices  in  the  closed  ( A  and  C)  and  modeled  20°-opened  ( B  and  D )  KcsA  structures  (M2  helices,  dark  gray.  Ml 
helices,  light  gray).  The  P  regions  of  the  channel  are  omitted  for  clarity.  Views  in  A  and  B  are  from  the  extracellular  end  of  the  channel.  Views  in  C  and 
D  are  side-on  from  the  membrane  plane  such  that  the  interhelical  gaps  appeared  to  be  maximized.  ( E)  Channel  pore  diameter  given  by  the  protein  van  der 
Waals  surface  versus  position  along  the  channel  axis,  as  measured  by  the  program  HOLE  (Smart  et  al.,  1993),  for  several  M2  helix  rotations  (x-ray  model, 
0°).  The  protein  extends  from  0-60  A  with  the  extracellular  end  oriented  at  left.  The  narrow  region  from  ~  10-25  A  contains  the  selectivity  filter,  and  the 
naiTow  region  from  ~30-50  A  in  the  x-ray  structure  is  lined  with  aliphatic  side  chains  from  the  four  M2  helices.  ( F)  Superimposed  ribbon  representations 
of  the  channel  in  the  closed  (gray)  and  20°-open  conformation  (dark  gray)  as  viewed  from  the  membrane  plane  (spheres,  light  gray)  for  electrostatics 
calculations.  The  remaining  sections  corresponding  to  the  selectivity  filter,  extended,  and  turret  regions  (black)  are  common  to  both  structures.  Orientations 
in  panels  C,  D,  and  F  are  with  the  extracellular  end  at  top. 
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respect  to  the  coordinate  axes  for  all  channels  studied.)  The  pKa  values 
were  computed  using  the  linear  Poisson-Boltzmann  electrostatics  module 
within  UHBD.  Test  calculations  with  the  nonlinear  Poisson-Boltz¬ 
mann  method  revealed  that  the  modifications  of  the  ionization  states  were 
not  large  enough  to  justify  their  greater  computational  demand.  The  ion¬ 
ization  states  were  not  substantially  affected  by  either  the  presence  of  a 
single  on-axis  K+  ion  or  by  ionic  strength  in  the  range  0.15  to  1.50  M. 
Symmetric  pH  conditions  were  used,  and  all  calculations  were  earned  out 
at  298  K. 

For  computing  potential  profiles  of  a  single  ion  along  the  channel  axis, 
the  calculated  ionization  states  were  converted  into  new,  effective  partial 
charges  for  the  side  chain  atoms.  The  total  force  acting  on  the  ion  as  a 
function  of  position  along  the  axis  was  then  calculated  using  UHBD.  The 
axial  component  of  this  force  was  integrated  over  the  length  of  the  channel 
to  give  the  potential  profile  of  the  ion.  Because  of  the  asymmetry  of  charge 
in  the  channel  protein  combined  with  the  absence  of  shielding  charges  in 
the  electrolyte  in  this  reduced  system,  it  was  necessary  to  normalize  the 
computed  transmembrane  potential  to  zero  by  the  addition  of  a  constant 
field. 


Molecular  dynamics:  determining  the  mobility  of 
ions  within  the  channel 

MD  simulations  of  ion  motions  throughout  the  channel  were  performed 
using  the  program  GROMOS96  (van  Gunsteren  et  al.,  1996)  with  solvation 
of  the  refined  version  of  the  channel  by  methods  we  previously  used  for  the 
gramicidin  channel  (Chiu  et  al.,  1993)  and  for  part  of  the  permeation 
pathway  for  a  model  of  the  sodium  channel  (Singh  et  al.,  1996).  Specifi¬ 
cally,  the  wide  parts  of  the  channel  pore  were  hydrated,  and  the  ends  of  the 
channel  were  bathed  in  approximately  hemispherical  caps  of  water  extend¬ 
ing  to  the  aromatic  collars.  Three  potassium  ions,  two  of  which  were 
initially  located  in  the  outer  crystallographic  positions  at  either  end  of  the 
selectivity  filter  and  one  ion  located  in  the  central  cavity,  were  included  for 
equilibration  for  1  ns  at  298  K.  To  prepare  the  system  for  ion  diffusion 
measurements,  two  of  the  ions  were  replaced  with  water  molecules,  and  the 
system  was  re-equilibrated  while  strongly  restraining  the  restraining  ion  to 
the  channel  axis.  The  position  of  the  strong  restraints  was  updated  in  1-A 
increments  along  the  channel  axis,  allowing  for  reequilibration  in  each 
instance.  Several  reequilibrated  systems  were  selected  to  represent  a  subset 
of  ion  positions  along  the  channel  axis,  including  those  in  which  the  ion 
was  located  in  the  outer  vestibule,  within  the  selectivity  filter,  in  the  central 
cavity,  and  in  the  narrow,  hydrophobic  part  of  the  channel.  Production  runs 
of  27  ps  with  no  restraints  on  the  ion  were  generated  in  each  region  of  the 
channel. 

Velocity  autocorrelation  functions  were  calculated  for  the  unrestrained 
ion  trajectories  in  each  part  of  the  channel.  The  velocity  autocorrelation 
function  (McQuarrie,  1976)  relevant  to  one-dimensional  motion  of  ions 
along  the  channel  axis  (see  below)  is  given  by  the  expression  (vz(?1)vz(?2)), 
in  which  vz  is  the  z-component  of  the  velocity  of  the  ion  at  two  times,  tx  and 
t2,  and  the  angle  brackets  denote  an  equilibrium  average.  In  an  equilibrium 
system  the  correlations  depend  on  only  differences  in  time  t  =  t2  —  t{.  We 
fit  the  velocity  autocorrelation  function  to  the  single  exponential  function 

<vz(0)vz(f)>  =  A0[v(0)]2exp(-f/r),  (1) 

in  which  A0  is  a  scaling  factor  and  the  time  constant,  t,  is  related  to  the  fluid 
dynamic  diffusion  coefficient,  D{,  by  the  expression 

r  =  mDf/kBT  (2) 

in  which  m  is  the  mass  of  the  ion,  T  is  the  temperature,  and  kB  is 
Boltzmann’s  constant. 
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TABLE  1  Input  parameters  for  the  Brownian  dynamics 
simulation  program  for  the  fully  opened  channel 


Parameter 

Value 

Channel  length 

6.0  nm 

Channel  mouth  diameter  (open  channel) 

1.35  nm 

Dielectric  constant  for  ion-ion  interactions 

20 

Channel  diffusion  coefficient 

2.51  X  10  5  cm2  s—  1 

Bath  diffusion  coefficient 

2.51  X  10-5  cm2 

Size  of  time  step 

0.4  ps 

Length  of  each  simulation 

0.2  ms 

Brownian  dynamics:  calculating  fluxes  and 
distributions  within  the  channel 

BD  simulations  were  performed  using  methods  we  have  previously  re¬ 
ported  (Cooper  et  al.,  1985;  Bek  and  Jakobsson,  1994)  to  compute  ion 
fluxes.  These  calculations  are  one-dimensional,  in  effect  constraining  the 
ions  to  move  along  the  channel  axis.  For  each  time  step,  each  ion  in  the 
system  is  moved  a  certain  distance  based  on  a  combination  of  deterministic 
motion  due  to  a  free  energy  gradient  and  random  motion  representing 
thermal  fluctuations,  according  to  the  equation 

D  (dU\ 

Az=  +  (3) 

in  which  Az  is  the  distance  that  an  ion  moves  in  the  time  interval  A t,  D  is 
the  ion  diffusion  coefficient,  R  is  the  ideal  gas  constant,  T  is  the  absolute 
temperature,  dU/dz  is  the  free  energy  gradient,  and  £  is  a  random  number 
chosen  from  a  normal  Gaussian  distribution  centered  at  zero.  The  free 
energy  gradient,  dU/dz,  is  approximated  by  its  electrostatic  component 
ZioirF(dV/dz),  in  which  ZiDn  is  the  valence  of  the  ion,  F  is  Faraday’s  constant, 
and  dV/dz  is  the  voltage  gradient.  The  voltage  gradient  is  the  force  seen  by 
an  ion  due  to  the  channel  and  surrounding  solvent  as  it  moves  through  the 
channel,  plus  the  applied  transmembrane  voltage,  plus  the  electric  field  due 
to  other  ions  in  the  channel. 

A  summary  of  the  parameters  for  the  BD  calculations  is  listed  in 
Table  1.  Based  on  the  MD  results  (described  in  the  Results  section 
below),  we  set  the  K+  diffusion  coefficient  inside  and  outside  the 
channel  equal  to  the  dynamic  ion  diffusion  coefficient  for  K+  in  bulk 
water.  A  dielectric  constant  of  20  was  used  for  ion-ion  interactions, 
which  is  consistent  with  the  ions  being  in  a  narrow  cavity  surrounded 
by  protein.  The  time  step  was  0.4  ps,  and  the  time  step  and  mobility 
together  lead  to  a  mean  thermal  jump  distance  of  0.45  A  per  time  step. 
Ions  enter  the  channel  by  means  of  the  “entrance  tube”  algorithm 
(Cooper  et  al.,  1985;  Jakobsson  and  Chiu,  1987;  Bek  and  Jakobsson, 
1994),  a  particular  implementation  of  the  flux-over-population  method 
(Farkas,  1927,  as  cited  by  Hanggi  et  al.,  1990).  The  capture  radius  for 
ions  attempting  to  enter  the  channel  was  taken  as  the  average  radius  of 
the  intracellular  half  of  the  channel  (Fig.  IE).  Thus  the  BD  simulations 
contain  no  arbitrarily  chosen  parameters. 

RESULTS  AND  DISCUSSION 
Computing  potential  profiles 

An  overview  of  our  open  channel  models  is  shown  in  Fig. 
1,  A-F.  Detailed  information  about  the  generation  of  these 
structures  is  given  in  the  Materials  and  Methods  section. 
Given  these  models,  we  first  compute  the  ionization  states 
of  the  titratable  residues  in  both  the  x-ray  and  opened 
channel  structures  at  pH  4  and  7,  as  it  has  been  demon- 
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TABLE  2  Comparison  of  ionization  states  of  titratable 
side  chains 


Residue 

pH  7 

pH  4 

GROMOS96 

0° 

20° 

0° 

20° 

His-25 

0.043 

0.021 

0.922 

0.931 

0 

Arg-27 

1.000 

1.000 

1.000 

1.000 

1 

Tyr-45 

0.000 

-0.001 

0.000 

0.000 

0 

Glu-51 

-0.998 

-0.997 

-0.598 

-0.486 

-1 

Arg-52 

1.000 

1.000 

1.000 

1.000 

1 

Tyr-62 

-0.003 

-0.003 

0.000 

0.000 

0 

Arg-64 

1.000 

1.000 

1.000 

1.000 

1 

Glu-7 1 

-0.486 

-0.513 

-0.017 

-0.020 

-1 

Tyr-78 

0.000 

0.000 

0.000 

0.000 

0 

Asp-80 

-0.999 

-0.999 

-0.973 

-0.974 

-1 

Tyr-82 

0.000 

0.000 

0.000 

0.000 

0 

Arg-89 

1.000 

1.000 

1.000 

1.000 

1 

Cys-90 

-0.001 

-0.002 

0.000 

0.000 

0 

Arg-117 

1.000 

1.000 

1.000 

1.000 

1 

Glu-118 

-0.963 

-0.997 

-0.315 

-0.354 

-1 

The  values  (in  units  of  e)  were  averaged  over  the  four  subunits.  The  data 
are  for  the  x-ray  (0°)  and  opened  channel  structures  (20°)  embedded  in  a 
dielectric  slab  and  exposed  to  electrolyte  (ionic  strength,  200  mM).  Also 
shown  are  the  standard  partial  charges  used  in  the  molecular  dynamics 
simulations  for  the  x-ray  structure  at  neutral  pH. 
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strated  experimentally  that  the  channel  is  functionally 
closed  at  pH  7  and  is  fully  open  at  pH  4  (Cuello  et  al.,  1998; 
Heginbotham  et  al.,  1999).  The  results  in  Table  2  show  that 
the  glutamates,  located  in  the  vestibule  and  toward  the  C 
terminus,  undergo  partial  protonation  at  low  pH,  whereas 
the  histidines,  toward  the  N  terminus,  carry  a  full  charge. 
(For  comparison,  the  values  used  by  the  GROMOS96  sim¬ 
ulation  package  standard  force  field  are  also  given  and  are 
in  close  correspondence  with  the  computed  ionization  states 
at  neutral  pH  for  all  residues  except  for  Glu-7 1 .)  Ionization 
states  of  side  chains  at  various  ionic  strengths  were  com¬ 
puted  similarly.  We  found  that  neither  the  channel  confor¬ 
mation  nor  variation  of  the  ionic  strength  had  an  appreciable 
effect  on  computed  ionization  states. 

Fig.  2  shows  the  potential  profiles  as  seen  by  a  K+  ion 
traversing  the  channel  along  the  channels  axis.  The  potential 
profile  for  the  x-ray  conformation  at  pH  7  has  three  notable 
features:  a  deep  potential  well  containing  the  selectivity 
filter,  a  second  potential  well  near  the  intracellular  end  of 
the  channel,  and  a  pronounced  potential  barrier  in  the  nar¬ 
row  region  of  the  channel.  All  curves  show  a  deep  potential 
well  near  the  extracellular  end  of  the  channel  that  is  due  to 
the  negatively  charged  carbonyl  oxygens  in  the  selectivity 
filter  and,  to  a  lesser  extent,  to  the  negatively  charged 
residues  in  the  vestibule  of  the  channel.  This  well  is  some¬ 
what  shallower  at  pH  4  than  at  pH  7  due  to  the  partial 
neutralization  of  Glu-51  in  the  vestibule.  The  potential  well 
near  the  intracellular  end  of  the  x-ray  structure  is  due  to 
negatively  charged  Glu-118  and  shows  a  large  pH  depen¬ 
dence.  The  pH  dependence  here  for  the  open  structure  is 
much  less  because  Glu-118  is  now  farther  from  the  channel 


FIGURE  2  (A)  Potential  profile  of  a  potassium  ion  in  translocating  the 
channel  along  the  channel’s  axis  in  the  x-ray  and  open  channel  structures 
at  neutral  and  acidic  pH  using  partial  charge  data  from  Table  2.  Those 
profiles  in  bold  relief  correspond  to  physiologically  relevant  conditions  and 
are  the  basis  for  the  subsequent  calculations.  ( B )  Potential  profiles  for  K+ 
in  the  structures  of  intermediate  open  sizes  (as  measured  by  the  degree  of 
M2  rotation)  at  pH  4.  Note  the  decline  in  the  Born  energy  barrier  near  the 
intracellular  end  of  the  channel  as  the  channel  is  opened. 

axis.  The  barrier,  resulting  from  the  dielectric  contrast  be¬ 
tween  the  high-dielectric  electrolyte  and  the  low-dielectric 
protein  and  surrounding  membrane,  is  an  image  potential 
barrier  (Sancho  et  al.,  1995)  that  can  be  thought  of  as  a  Born 
energy,  or  desolvation  energy,  associated  with  moving  an 
ion  from  the  surrounding  electrolyte  into  the  narrow  non¬ 
polar  channel  interior.  The  profile  for  the  opened  channel  at 
pH  4  reveals  that  this  image  potential  barrier  is  greatly 
reduced  upon  widening  the  channel  lumen.  It  is  seen  that  the 
depth  of  the  potential  well  at  the  selectivity  filter  is  not  very 
sensitive  to  the  degree  to  which  the  channel  is  opened  but  is 
sensitive  to  the  ionization  state  of  the  protein.  The  potential 
profiles  as  computed  for  other  ionic  strengths  (not  shown) 
provide  similar  results.  Fig.  2  B  shows  the  potential  profile 
for  the  potassium  ion  traversing  the  channel  axis  for  differ¬ 
ent  sizes  of  channel  openings  at  pH  4.  It  is  seen  that  the 
potential  barrier  near  the  intracellular  end  of  the  channel  is 
systematically  reduced  as  the  channel  is  opened  further. 

In  summary,  our  modeling  and  electrostatics  calculations 
thus  show  that  the  depth  of  the  potential  well  at  the  selec¬ 
tivity  filter  is  controlled  by  the  ionization  states  of  titratable 
residues  of  the  protein,  whereas  the  image  potential  barrier 
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near  the  intracellular  half  of  the  channel  is  controlled  by  the 
size  of  the  pore  at  the  intracellular  end. 


Computing  ion  mobilities 

Velocity  autocorrelation  functions  were  calculated  from  the 
ion  trajectories  generated  by  MD  simulations.  Values  for  the 
fluid  dynamic  diffusion  coefficient  (due  to  local  friction), 
determined  from  Eqs.  1  and  2,  were  found  to  range  from 
1.8  X  1CT5  to  3.3  X  10~5  cm2  s-1,  in  close  agreement  with 
the  experimental  value  of  2.5  X  1 0  5  cm2  s  !  that  pertains 
to  K+  ions  in  bulk  water  (Hille,  1992).  Because  of  the 
complexity  of  the  ion  solvation  environment,  the  velocity 
autocorrelation  functions  were  more  complex  than  simple 
exponentials.  However,  the  best  fit  to  an  exponential  still 
suffices  to  give  an  approximately  correct  diffusion  coeffi¬ 
cient.  In  the  wide  part  of  the  channel,  the  ion  is  in  fact 
solvated  by  water.  In  the  narrow  part  of  the  channel,  it 
appears  from  our  results  that  the  carbonyl  oxygens  substi¬ 
tute  approximately  equivalently  to  bulk  water  oxygens  in 
determining  the  local  friction  for  ion  movement.  It  should 
be  pointed  out  that  the  diffusion  coefficient  in  this  case  is 
not  the  effective  diffusion  coefficient,  but  rather  the  recip¬ 
rocal  of  the  local  friction.  The  effective  diffusion  coefficient 
in  the  selectivity  filter  is  certainly  less  than  in  bulk  water, 
due  to  the  potential  well  created  by  the  carbonyl  oxygens, 
even  while  the  local  friction  is  the  same  as  bulk  water.  This 
situation  is  similar  to  that  pertaining  in  gramicidin,  where 
we  earlier  found  that  the  local  friction  for  sodium  ions  in 
gramicidin  was  similar  to  bulk  water,  whereas  the  effective 
diffusion  coefficient  was  lower  by  a  factor  of  10,  due  to  the 
potential  wells  created  by  carbonyl  oxygens  in  that  structure 
(Chiu  et  al.,  1993). 

Computing  ionic  fluxes 

Ionic  fluxes  were  calculated  using  a  one-dimensional  BD 
method  for  the  ions.  The  resulting  conductance  calculated 
for  the  closed  channel  at  both  pH  7  and  4  was  a  fraction  of 
a  picosiemen.  In  Fig.  3  A  the  computed  current-voltage 
(I-V)  curve  for  the  fully  opened  channel  for  250  mM  K+ 
shows  an  approximately  linear  relationship  over  the  voltage 
range  of  ±100  mV  in  agreement  with  the  experiment.  The 
corresponding  average  channel  conductance  over  this  volt¬ 
age  range  is  — 110  pS,  in  reasonable  agreement  with  exper¬ 
imentally  determined  maximal  conductances  of  90  pS 
(Schrempf  et  al.,  1995)  and  135  pS  at  200  mM  K+  (Cuello 
et  al.,  1998)  and  of  —130  pS  at  250  mM  K+  (Meuser  et  al., 
1999).  The  mean  number  of  ions  in  the  channel  and  the 
Ussing  flux  ratio  exponent  were  computed  for  various  trans¬ 
membrane  voltages  (Table  3)  and  were  found  to  be  close  in 
value  to  each  other,  in  agreement  with  the  original  theory 
relating  these  two  quantities  (Hodgkin  and  Keynes,  1955). 
Although  there  are  currently  no  experimental  flux  ratio 
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FIGURE  3  (A)  Current-voltage  curves  as  computed  by  BD  for  the 

opened  channel  at  pH  4  for  250  mM  K+  using  the  parameters  in  Table  1, 
compared  with  experimental  results  (Meuser  et  al.,  1999)  for  the  largest 
conducting  substate  under  comparable  conditions.  Each  simulated  point 
represents  a  0.2-ms-long  simulation.  (B)  Current-concentration  curves  for 
the  opened  state  at  pH  4  for  several  membrane  potentials  with  parameters 
from  Table  1,  compared  with  experiment  (Meuser  et  al.,  1999).  Each  curve 
is  normalized  to  the  respective  value  of  the  current  at  250  mM  K+.  (C) 
Current- voltage  curves  as  computed  by  BD  for  channels  of  different 
opened  sizes  at  250  mM  K+.  Two  factors  are  responsible  for  the  difference 
in  currents.  One  is  the  size  of  the  Born  potential  barrier  near  the  intracel¬ 
lular  end  of  the  channel,  as  shown  in  Fig.  2  B.  The  other  is  the  capture 
diameter  for  ions  to  impinge  on  the  intracellular  end  of  the  channel,  as 
visualized  in  Fig.  1  E.  The  capture  diameters  used  for  the  calculations 
shown  were  5.5  A  for  5°,  8.4  A  for  10°,  1 1.0  A  for  15°,  and  13.5  A  for  20°. 
In  all  calculations  shown  in  Fig.  3  the  standard  deviations  of  the  simulated 
currents  are  smaller  than  the  size  of  the  data  symbols  used. 
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TABLE  3  Mean  number  of  ions  in  the  channel  and  Ussing 
flux  ratio  exponent  as  a  function  of  transmembrane  voltage 
for  the  opened  channel  at  pH  4  and  at  250  mM  K+,  as 
computed  by  Brownian  dynamics 


Voltage  (mV) 

Mean  number 

Ussing  flux  ratio  exponent 

100 

2.01 

2.08 

80 

2.01 

2.22 

60 

2.02 

2.16 

40 

2.02 

2.25 

20 

2.03 

2.14 

0 

2.02 

Indeterminate 

-20 

2.03 

2.00 

-40 

2.03 

2.16 

-60 

2.03 

2.21 

-80 

2.02 

2.09 

-100 

2.01 

2.01 

exponent  data  for  KcsA,  a  recent  study  (Stampe  et  al.,  1998) 
on  an  inward-rectifying  K  channel  has  shown  ratios  be¬ 
tween  2.1  and  2.5  for  membrane  potentials  ranging  from 
—50  to  —25  mV,  which  is  in  close  agreement  with  our 
simulations. 

Fig.  3  B  shows  the  relationship  between  current  (normal¬ 
ized  to  the  250  mM  K+  value)  and  concentration  for  the 
opened  channel  at  several  transmembrane  voltages  with 
comparison  with  experimental  values  (Meuser  et  al.,  1999). 
(The  channel-ion  potential  profiles  were  recalculated  for 
each  ion  concentration.)  Whereas  the  calculated  concentra¬ 
tion  dependence  of  the  simulated  current  is  similar  to  that 
seen  in  experiment,  the  simulated  currents  tend  to  saturate  at 
a  slightly  lower  concentrations.  We  believe  this  behavior 
may  be  due  to  the  approximation  of  one-dimensional  ion 
motion  implicitly  used  in  the  BD  simulations,  an  approxi¬ 
mation  that  is  not  as  good  at  near-saturating  concentrations 
as  it  is  at  lower  more  physiological  concentrations. 

Fig.  3  C  shows  the  I-V  curves  for  single  channels  for  the 
different  open  sizes  shown  in  Fig.  1  E.  The  modeled  par¬ 
tially  opened  channels  produce  intermediate  conductances 
in  a  fashion  similar  to  conducting  substates.  The  conduc¬ 
tance  increase  as  the  pore  size  increases  is  partly  due  to  the 
reduction  in  the  Born  potential  barrier  near  the  intracellular 
end  (Fig.  2  B )  and  partly  due  to  the  increase  in  the  capture 
diameter  of  the  channel  when  the  channel  opening  is  larger. 
It  may  be  that  the  basis  of  conducting  substates  in  K 
channels  is  the  size  of  the  pore  aperture  at  the  intracellular 
end  of  the  channel  in  a  similar  fashion  to  this  model. 

Features  of  Brownian  ion  trajectories 

The  sample  of  BD  ion  trajectories  in  Fig.  4  demonstrates  the 
extreme  ability  of  the  negative  charges  lining  the  selectivity 
filter  to  concentrate  cations.  Here  the  ribbon  structure  of  the 
opened  channel  is  superposed  on  the  trajectories  so  that  the 
structural  basis  for  the  preferred  ionic  positions  can  be  seen. 
There  is  a  clear  preference  for  a  double-occupancy  state  in 


the  selectivity  filter,  consistent  with  the  computed  mean 
number  of  ions  in  the  channel  (Table  3).  The  histogram  of 
ion  positions  computed  for  the  entire  length  of  the  BD 
simulation  reveals  two  prominent  peaks  in  the  selectivity 
filter  separated  by  ~7  A,  similar  to  the  pattern  of  ion 
distribution  as  revealed  by  x-ray  crystallography  (Doyle  et 
al.,  1998).  Fig.  4  also  shows  the  occasional  appearance  of  a 
single-occupancy  state  in  the  selectivity  filter  where  the 
mean  position  of  the  single  ion  lies  in  the  carbonyl  “cage” 
in  between  the  two  that  are  preferred  during  double  occu¬ 
pancy.  This  relationship  between  singly  and  doubly  occu¬ 
pied  states  is  similar  to  the  findings  of  recent  free-energy 
perturbation  MD  simulations  (Aqvist  and  Luzhkov,  2000). 

The  trajectories  in  Fig.  4  also  reveal  that  two  ions  occu¬ 
pying  the  selectivity  filter  move  nearly  in  unison  in  response 
to  an  incoming  third  ion.  The  movement  of  one  ion  out  of 
one  end  of  the  selectivity  filter  and  the  entrance  of  another 
from  the  other  end  is  so  close  to  simultaneous  that  on  the 
time  scale  of  Fig.  4,  the  trajectories  suggest  a  “knock-off' 
mechanism.  Recent  experimental  data  on  channel  blocking 
of  KcsA  by  barium  ions  has  suggested  similarly  that  a  third 
ion  could  enter  the  selectivity  filter  to  push  the  queue  of  ions 
along  (Jiang  and  MacKinnon,  2000). 

Because  the  selectivity  filter  holds  two  ions  practically  all 
the  time,  and  because  the  mean  channel  occupancy  is  only 
fractionally  over  two  (Table  3),  it  follows  that  the  wide  part 
of  the  channel  hardly  ever  holds  more  than  one  ion.  This 
shows  why  the  one-dimensional  approximation  to  the  mo¬ 
tion  is  reasonably  accurate.  In  the  selectivity  filter  the  mo¬ 
tion  really  is  one-dimensional  because  the  channel  is  narrow 
there.  In  the  wide  part  of  the  channel,  only  the  projection  of 
the  motion  of  the  single  ion  onto  the  channel  axis  is  all  that 
contributes  to  the  flux.  Thus,  ignoring  the  lateral  motion  of 
the  ions  in  this  region  does  not  introduce  gross  inaccuracies 
into  the  calculation. 

Ion  channels  are  essentially  electrical  resistance  elements 
(Hodgkin  and  Huxley,  1952).  In  the  KcsA  (and  presumably 
other)  potassium  channels,  the  selectivity  filter  is  in  series 
with  the  wider  portions  of  the  channel  at  the  extracellular 
and  intracellular  ends.  We  analyzed  this  channel  as  two 
resistors  in  series,  one  comprising  the  selectivity  filter  and 
the  other  the  rest  of  the  channel.  Because  the  selectivity 
filter  is  almost  always  doubly  occupied,  regardless  of  bath 
electrolyte  concentration,  we  postulate  that  the  resistance  of 
the  selectivity  filter  is  a  constant.  We  further  postulate  that 
the  resistance  of  the  rest  of  the  channel  is  inversely  propor¬ 
tional  to  electrolyte  concentration,  similar  to  bulk  electro¬ 
lyte.  This  leads  to  the  functional  form 

R  =  Rf  +  A/[  K+],  (4) 

in  which  R  is  the  resistance  of  the  single  channel  (reciprocal 
of  the  conductance),  Rt  is  the  intrinsic  resistance  of  the 
selectivity  filter,  and  A  is  a  constant.  This  model  of  resistors 
in  series  provides  a  physical  basis  for  the  frequently  ob- 
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FIGURE  4  A  120-ns  sample  of  BD  ion  trajectories  in  the  opened  channel  at  pH  4  at  0  mV  membrane  potential  in  250  mM  K+  solution  using  parameters 
in  Table  1.  For  ease  of  visualization  the  data  were  smoothed  using  a  window  of  120  ps,  effectively  filtering  out  oscillations  greater  than  ~1010  Hz.  Ions 
that  have  successfully  traversed  the  channel  in  this  time  window  are  shown  in  bold  relief.  A  RasMol  (Sayle  and  Milner- White,  1995)  ribbon  representation 
of  the  backbones  (gray)  from  two  opposing  subunits  showing  the  negatively  charged  carbonyl  groups  (black)  projecting  into  the  selectivity  filter  lumen 
is  superposed. 


served  property  of  ion  channels  to  follow  Michaelis-Men- 
ten  kinetics  (Aidley  and  Stanfield,  1996).  The  solid  and 
dashed  lines  in  Fig.  5  show  the  fit  of  Eq.  4  to  both  computed 
and  experimental  data  (same  data  as  Fig.  3  B )  for  a  voltage 
of  100  mV.  It  is  seen  that  the  fit  is  very  good.  Similarly 
good  fits  to  Eq.  4  are  obtained  for  the  full  range  of  simulated 
voltages.  In  this  view,  the  reciprocal  of  Rt  (—500  pS)  is  the 
theoretical  maximal  conductance  achievable  by  a  potassium 
channel  if  all  resistance  in  series  with  the  selectivity  filter 
could  be  eliminated.  The  ratio  of  selectivity  filter  resistance 
to  that  of  the  wider  pore  regions  was  computed  (Fig.  5, 
inset).  At  low  bath  potassium  concentrations,  the  majority 
of  the  KcsA  resistance  appears  to  be  in  the  wider  regions  of 
the  channel,  whereas  at  high  concentrations  the  majority  of 
the  resistance  is  in  the  selectivity  filter. 

It  may  at  first  seem  counterintuitive  that  at  physiologi¬ 
cally  relevant  ion  concentrations  the  majority  of  the  resis¬ 
tance  is  in  the  wide  part  of  the  channel.  To  understand  this, 
consider  that  the  conductance  of  any  medium  through  which 
ions  flow  is  a  product  of  the  mobility  and  the  charge  density. 
To  some  extent  the  effective  mobility  of  the  ions  is  reduced 
in  the  selectivity  filter  by  the  negative  charges  of  the  back¬ 
bone  carbonyl  oxygens,  yet  the  concentration  is  increased 


FIGURE  5  Opened-channel  resistance,  R,  computed  from  the  ratio  of 
voltage  to  current,  as  a  function  of  electrolyte  concentration  at  100  mV  for 
BD  simulations  (circles)  and  for  electrophysiological  measurements  of  the 
maximal  conducting  substate  (Fig.  4  of  Meuser  et  al.,  1999)  (triangles). 
Fitting  of  the  data  sets  to  Eq.  4  yielded  constants  of  R{  =  2.3  Gfl  and  A  = 
1.34  X  109  fFM  for  the  simulations  and  Rf  =  1.6  Gfl  and  A  =  1.48  X  109 
fl*M  for  the  experiments.  The  inset  shows  the  ratio  of  the  resistance  of  the 
wide  regions  of  the  channel  to  the  resistance  of  the  selectivity  filter  for  the 
simulations. 
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enormously  (~20  M  K+),  which  is  more  than  compensating 
for  the  reduction  in  effective  mobility. 

SUMMARY  AND  CONCLUSIONS 

We  have  described  and  carried  out  a  hierarchical  strategy 
for  characterizing  the  permeation  of  the  KcsA  potassium 
channel.  This  strategy  begins  with  the  x-ray  crystal  structure 
of  the  channel  (Doyle  et  al.,  1998)  from  which  a  series  of 
opened  channel  structures  was  built  by  a  process  of  rotating 
the  transmembrane  helices  surrounding  the  channel  perme¬ 
ation  pathway  in  a  manner  generally  consistent  with  block¬ 
ing  experiments,  toxin-binding  studies,  and  hydrophobic 
matching  principles.  Experimental  pH  conditions  for  the 
functionally  closed  and  open  conformations  were  incorpo¬ 
rated  through  pKa  calculations  of  the  ionization  states  of  the 
titratable  residues.  Electrostatics  calculations  were  used  to 
obtain  the  potential  profiles  of  ions  along  the  channel  axis. 
MD  simulations  were  performed  to  obtain  the  ion  diffusion 
coefficient  for  incorporation  into  BD  calculations  for  gen¬ 
erating  ion  fluxes.  Although  the  calculations  have  many 
approximations  (e.g.,  representing  the  ion-channel  potential 
profile  by  an  electrostatic  calculation,  representing  the  ion 
motion  through  the  channel  as  one-dimensional,  represent¬ 
ing  ion-ion  interaction  in  the  channel  as  governed  by  an 
equivalent  dielectric  constant,  not  having  a  high-resolution 
structure  for  the  open  form  of  the  channel),  the  calculations 
appear  to  capture  the  essence  of  the  permeation  through  the 
KcsA  channel,  based  on  the  close  correspondence  between 
the  behavior  of  the  real  channel  and  the  BD  simulations. 

Examination  of  the  calculations  suggests  several  features 
of  K+  permeation  through  potassium  channels.  1)  The  x-ray 
model  structure  of  the  KcsA  channel  is  closed  not  by 
physical  occlusion  but  rather  by  image  forces  that  arise  in 
attempting  to  move  an  ion  from  the  electrolyte  to  a  narrow, 
nonpolar  cavity.  2)  Changing  the  ionization  states  of  the 
titratable  residues  by  reducing  the  pH  from  7  to  4  in  the 
x-ray  crystal  structure  does  not  alter  the  channel-ion  poten¬ 
tial  profde  sufficiently  to  account  for  the  increase  in  con¬ 
ductance  in  lowering  the  pH.  Therefore  a  conformational 
change  in  the  transmembrane  helices  of  the  channel  is 
required  for  ion  conduction.  3)  The  local  friction  of  ions 
moving  within  the  channel  is  similar  to  that  of  ions  in  bulk 
water.  However,  the  effective  mobility  in  the  selectivity 
filter  is  smaller  due  to  the  partial  confinement  of  ions  in  the 
selectivity  filter  by  negative  structural  charges  in  the  protein 
(Chiu  et  al.,  1993;  Allen  et  al.,  1999).  This  can  be  seen  in  the 
computed  trajectories  in  Fig.  4,  noting  that  the  fluid  dy¬ 
namic  diffusion  coefficient  for  ions  in  the  selectivity  filter  is 
the  same  as  for  ions  elsewhere  in  the  channel,  but  the  ions 
in  the  selectivity  filter  are  effectively  immobilized  by  the 
negative  charges  in  the  carbonyl  oxygens  lining  the  selec¬ 
tivity  filter  lumen.  4)  The  most  common  occupancy  state  of 
the  open  channel  at  near-physiological  concentrations  is  two 
ions  in  the  selectivity  filter.  In  and  near  the  selectivity  filter 


are  where  ion-ion  interactions  primarily  occur.  The  domi¬ 
nant  mode  for  ion  motion  through  the  selectivity  filter  is 
effectively  by  a  knock-off  mechanism  involving  three  ions, 
as  suggested  experimentally  (Jiang  and  MacKinnon,  2000). 
5)  For  near-physiological  electrolyte  concentrations  there  is 
seldom  multiple  occupancy  in  the  extracellular  vestibule  or 
along  the  intracellular  half  of  the  channel.  Ion  motion  in 
those  regions  is  therefore  largely  governed  by  the  electro¬ 
diffusion  of  single  ions.  It  is  for  this  reason  that  the  one¬ 
dimensional  approximation  to  the  motion  embodied  in  the 
Brownian  computational  model  in  this  paper  gives  reason¬ 
ably  accurate  estimates  of  the  flux.  6)  The  overall  saturation 
behavior  of  the  channel  in  symmetric  solutions  is  rather  well 
fit  by  a  simple  model  of  two  resistors  in  series,  one  for  the 
selectivity  filter  and  the  other  for  the  wider  parts  of  the 
channel.  The  resistance  of  the  selectivity  filter  is  ~2  X  109 
(1.  The  combined  resistance  of  the  wider  parts  of  the  chan¬ 
nel  is  approximately  inversely  proportional  to  the  bath  con¬ 
centration  of  permeant  ions. 

Our  calculations  are  generally  consistent  with  other  com¬ 
putational  studies  of  the  KcsA  channel.  Our  MD  trajectories 
for  ions  in  the  selectivity  filter  and  in  the  wider  regions  of 
the  channel  in  general  agree  with  those  of  other  workers 
(Berneche  and  Roux,  2000;  Guidoni  et  al.,  2000;  Shrivas- 
tava  and  Sansom,  2000)  The  observation  in  the  BD  calcu¬ 
lations  that  the  dominant  occupancy  state  in  the  channel  is 
two  ions  ~1  A  apart  agrees  with  free  energy  perturbation 
studies  (Aqvist  and  Luzhkov,  2000),  as  well  as  being  also 
suggested  by  the  x-ray  model  structure  (Doyle  et  al.,  1998). 
Our  ionization  state  for  the  residue  Glu-71  is  somewhat 
more  negatively  charged  than  suggested  by  others  (Roux 
and  MacKinnon,  1999;  Aqvist  and  Luzhkov,  2000),  perhaps 
due  to  a  somewhat  different  methodology  of  computing  the 
p Ka.  As  these  residues  are  over  6  A  from  the  channel  axis 
and  are  largely  occluded  by  the  selectivity  filter  from  the 
channel  axis,  this  difference  is  unlikely  to  affect  the  con¬ 
clusions  of  this  study.  We  also  appear  to  show  somewhat 
less  tendency  for  an  ion  to  occupy  the  wide  part  of  the 
channel  than  previously  suggested  (Roux  and  MacKinnon, 
1999).  This  difference  may  be  due  to  the  fact  that  our  ion 
distributions  are  computed  for  an  opened  channel.  Also, 
Roux  and  MacKinnon  based  their  main  conclusions  on 
calculations  in  which  the  protein  was  assigned  a  dielectric 
constant  of  2,  whereas  our  assumed  dielectric  constant  of  20 
would  somewhat  attenuate  the  effect  of  the  pore  helix 
potentials  to  attract  an  ion  into  the  cavity. 

In  the  approach  described  in  this  paper,  molecular  mod¬ 
eling  based  on  structural  data,  electrostatic  calculations, 
MD,  and  BD  simulations  are  combined  to  achieve  a  com¬ 
prehensive  view  of  the  permeation  process  in  the  KcsA  ion 
channel.  The  same  overall  approach  seems  capable  of  ap¬ 
plication  to  other  ion  channels  as  well. 

This  work  was  supported  by  the  National  Science  Foundation.  Computer 
time  from  the  National  Computational  Science  Alliance,  which  is  also 


37 


Biophysical  Journal  81(5)  2473-2483 


2482 


Mashl  et  al. 


largely  funded  by  the  National  Science  Foundation,  is  gratefully  acknowl¬ 
edged.  Simulations  were  earned  out  on  the  Origin  2000  at  the  National 
Center  for  Supercomputing  Applications  at  the  University  of  Illinois.  We 
thank  Dr.  Robin  Shealy  for  preparing  the  refined  KcsA  structure  and  Drs. 
Shankar  Subramaniam,  Larry  Scott,  and  Serdar  Kuyucak  for  helpful  com¬ 
ments. 


REFERENCES 

Aidley,  D.  J.,  and  P.  R.  Stanfield.  1996.  Ion  Channels:  Molecules  in  Action. 
Cambridge  University  Press,  Cambridge,  UK. 

Allen,  T.  W.,  A.  Bliznyuk,  A.  P.  Rendell,  S.  Kuyucak,  and  S.-H.  Chung. 
2000.  The  potassium  channel:  structure,  selectivity  and  diffusion. 
J.  Chem.  Phys.  112:8191-8204. 

Allen,  T.  W.,  S.  Kuyucak,  and  S.-H.  Chung.  1999.  Molecular  dynamics 
study  of  the  KcsA  potassium  channel.  Biophys.  J.  77:2502-2516. 

Antosiewicz,  J.,  J.  A.  McCammon,  and  M.  K.  Gilson.  1994.  Prediction  of 
pH-dependent  properties  of  proteins.  /.  Mol.  Biol.  238:415-436. 

Aqvist,  J.,  and  V.  Luzhkov.  2000.  Ion  permeation  mechanism  of  the 
potassium  channel.  Nature.  404:881-884. 

Bek,  S.,  and  E.  Jakobsson.  1994.  Brownian  dynamics  study  of  a  multiply- 
occupied  cation  channel:  application  to  understanding  permeation  in 
potassium  channels.  Biophys.  J.  66:1028-1038. 

Bemeche,  S.,  and  B.  Roux.  2000.  Molecular  dynamics  of  the  KcsA  K+ 
channel  in  a  bilayer  membrane.  Biophys.  J.  78:2900-2917. 

Brooks,  B.  R.,  R.  E.  Bruccoleri,  B.  D.  Olafson,  D.  J.  States,  S.  Swami- 
nathan,  and  M.  Karplus.  1983.  CHARMM:  a  program  for  macromolec- 
ular  energy,  minimization,  and  dynamics  calculations.  J.  Comp.  Chem. 
4:187-217. 

Chiu,  S.-W.,  J.  A.  Novotny,  and  E.  Jakobsson.  1993.  The  nature  of  ion  and 
water  barrier  crossings  in  a  simulated  ion  channel.  Biophys.  J.  64: 
98-109. 

Chung,  S.-H.,  T.  W.  Allen,  M.  Hoyles,  and  S.  Kuyucak.  1999.  Permeation 
of  ions  across  the  potassium  channel:  Brownian  dynamics  studies.  Bio¬ 
phys.  J.  77:2517-2533. 

Cooper,  K.,  E.  Jakobsson,  and  P.  Wolynes.  1985.  The  theory  of  ion 
transport  through  membrane  channels.  Prog.  Biophys.  Mol.  Biol.  46: 
51-96. 

Corry,  B.,  S.  Kuyucak,  and  S.-H.  Chung.  2000.  Tests  of  continuum  theories 
as  models  of  ion  channels:  II.  Poisson-Nernst-Planck  theory  versus 
Brownian  dynamics.  Biophys.  J.  78:2364-2381. 

Cowan,  S.  W.,  T.  Schirmer,  G.  Rummel,  M.  Steiert,  R.  Ghosh,  R.  A. 
Pauptit,  J.  N.  Jansonius,  and  J.  P.  Rosenbusch.  1992.  Crystal  structures 
explain  functional  properties  of  two  E.  coli  porins.  Nature.  358:727-733. 

Cuello,  L.  G.,  J.  G.  Romero,  D.  M.  Cortes,  and  E.  Perozo.  1998.  pH- 
dependent  gating  in  the  Streptomyces  lividans  K+  channel.  Biochemis¬ 
try.  37:3229-3236. 

Doyle,  D.  A.,  J.  M.  Cabral,  R.  A.  Pfuetzner,  A.  L.  Kuo,  J.  M.  Gulbis,  S.  L. 
Cohen,  B.  T.  Chait,  and  R.  MacKinnon.  1998.  The  structure  of  the 
potassium  channel:  molecular  basis  of  K+  conduction  and  selectivity. 
Science.  280:69-77. 

Farkas,  L.  1927.  Rate  of  nucleus  formation  in  saturated  vapors.  Z.  Physik. 
Chem.  (Leipzig).  125:236-242. 

Gibas,  C.  J.,  S.  Subramaniam,  J.  A.  McCammon,  B.  C.  Braden,  and  R.  J. 
Poljak.  1997.  pH  dependence  of  antibody /lysozyme  complexation.  Bio¬ 
chemistry.  36:15599-15614. 

Gilson,  M.  K.  1993.  Multiple-site  titration  and  molecular  modeling:  two 
rapid  methods  for  computing  energies  and  forces  for  ionizable  groups  in 
proteins.  Proteins.  15:266-282. 

Guidoni,  L.,  V.  Tone,  and  P.  Carloni.  2000.  Water  and  potassium  dynam¬ 
ics  inside  the  KcsA  K+  channel.  FEBS  Lett.  477:37-42. 

Hanggi,  P.,  P.  Talkner,  and  M.  Borkovec.  1990.  Reaction-rate  theory:  50 
years  after  Kramers.  Rev.  Mod.  Phys.  62:251-341. 


Heginbotham,  L.,  M.  LeMasurier,  L.  Kolmakova-Partensky,  and  C.  Miller. 
1999.  Single  Streptomyces  lividans  K+  channels:  functional  asymme¬ 
tries  and  sidedness  of  proton  activation.  J.  Gen.  Physiol.  114:551-559. 

Hille,  B.  1992.  Ionic  Channels  of  Excitable  Membranes,  2nd  ed.  Sinauer, 
Sunderland,  MA. 

Hodgkin,  A.  L.,  and  A.  F.  Huxley.  1952.  A  quantitative  description  of 
membrane  current  and  its  application  to  conduction  and  excitation  in 
nerve.  J.  Physiol.  117:500-544. 

Hodgkin,  A.  L.,  and  R.  D.  Keynes.  1955.  The  potassium  permeability  of  a 
giant  nerve  fibre.  J.  Physiol.  128:61-88. 

Huang,  C.-J.,  I.  Favre,  and  E.  Moczydlowski.  2000.  Permeation  of  large 
tetraalkyl  ammonium  cations  through  mutant  and  wild-type  voltage¬ 
gated  Na+  channels  as  revealed  by  relief  of  block  at  high  voltage.  J.  Gen. 
Physiol.  115:435-453. 

Jakobsson,  E.  1998.  Using  theory  and  simulation  to  understand  permeation 
and  selectivity  in  ion  channels.  Methods.  14:342-351. 

Jakobsson,  E.,  and  S.-W.  Chiu.  1987.  Stochastic  theory  of  ion  movement 
in  channels  with  single-ion  occupancy.  Biophys.  J.  52:33-45. 

Jiang,  Y.,  and  R.  MacKinnon.  2000.  The  barium  site  in  a  potassium 
channel  by  x-ray  crystallography.  J.  Gen.  Physiol.  115:269-272. 

Jorgensen,  W.  L.,  and  J.  Tirado-Rives.  1988.  The  OPLS  potential  functions 
for  proteins:  energy  minimizations  for  crystals  of  cyclic  peptides  and 
crambin.  J.  Am.  Chem.  Soc.  110:1657-1666. 

Killian,  J.  A.  1998.  Hydrophobic  mismatch  between  proteins  and  lipids  in 
membranes.  Biochim.  Biophys.  Acta.  1376:401-416. 

Laskowski,  R.  A.,  M.  W.  Mac  Arthur,  D.  S.  Moss,  and  J.  M.  Thornton. 
1993.  PROCHECK:  a  program  to  check  the  stereochemical  quality  of 
protein  structures.  J.  Appl.  Crystallogr.  26:283-291. 

Madura,  J.  D.,  J.  M.  Briggs,  R.  C.  Wade,  M.  E.  Davis,  B.  A.  Luty,  A.  Ilin, 
J.  Antosiewicz,  M.  K.  Gilson,  B.  Bagheri,  L.  R.  Scott,  and  J.  A. 
McCammon.  1995.  Electrostatics  and  diffusion  of  molecules  in  solution: 
simulations  with  the  University  of  Houston  Brownian  dynamics  pro¬ 
gram.  Comp.  Phys.  Commun.  91:57-95. 

McQuarrie,  D.  A.  1976.  Statistical  Mechanics.  Harper  Collins,  New  York. 
452-466. 

Meuser,  D.,  H.  Splitt,  R.  Wagner,  and  H.  Schrempf.  1999.  Exploring  the 
open  pore  of  the  potassium  channel  from  Streptomyces  lividans.  FEBS 
Lett.  462:447-452. 

Perozo,  E.,  D.  M.  Cortes,  and  L.  G.  Cuello.  1999.  Structural  rearrange¬ 
ments  underlying  K+  channel  activation  gating.  Science.  285:73-78. 

Roux,  B.,  and  R.  MacKinnon.  1999.  The  cavity  and  pore  helices  in  the 
KcsA  K+  channel:  electrostatic  stabilization  of  monovalent  cations. 
Science.  285:100-102. 

Sali,  A.,  J.  P.  Overington,  M.  S.  Johnson,  and  T.  L.  Blundell.  1990.  From 
comparisons  of  protein  sequences  and  structures  to  protein  modeling  and 
design.  Trends  Biochem.  Sci.  15:235-240. 

Sancho,  M.,  M.  B.  Partenskii,  V.  Dorman,  and  P.  C.  Jordan.  1995.  Ex¬ 
tended  dipolar  chain  model  for  ion  channels:  electrostriction  effects  and 
the  translocational  energy  barrier.  Biophys.  J.  68:427-433. 

Sayle,  R.  A.,  and  E.  J.  Milner- White.  1995.  RasMol:  biomolecular  graphics 
for  all.  Trends  Biochem.  Sci.  20:374-376. 

Schrempf,  H.,  O.  Schmidt,  R.  Kummerlen,  S.  Hinnah,  D.  Muller,  M. 
Betzler,  T.  Steinkamp,  and  R.  Wagner.  1995.  A  prokaryotic  potassium 
ion  channel  with  two  predicted  transmembrane  segments  from  Strepto¬ 
myces  lividans.  EMBO  J.  14:5170-5178. 

Shrivastava,  I.  H.,  C.  E.  Capener,  L.  R.  Forrest,  and  M.  S.  P.  Sansom.  2000. 
Structure  and  dynamics  of  K  channel  pore-lining  helices:  a  comparative 
simulation  study.  Biophys.  J.  78:79-92. 

Shrivastava,  I.  H.,  and  M.  S.  P.  Sansom.  2000.  Simulations  of  ion  perme¬ 
ation  through  a  potassium  channel:  molecular  dynamics  of  KcsA  in  a 
phospholipid  bilayer.  Biophys.  J.  78:557-570. 

Singh,  C.,  R.  Sankararamakrishnan,  S.  Subramaniam,  and  E.  Jakobsson. 
1996.  Solvation,  water  permeation,  and  ionic  selectivity  of  a  putative 
model  for  the  pore  region  of  the  voltage-gated  sodium  channel.  Biophys.  J. 
71:2276-2288. 

Smart,  O.  S.,  J.  M.  Goodfellow,  and  B.  A.  Wallace.  1993.  The  pore 
dimensions  of  gramicidin  A.  Biophys.  J.  65:2455-2460. 


Biophysical  Journal  81(5)  2473-2483 


38 


Predicting  Permeation  in  Ion  Channels 


2483 


Splitt,  H.,  D.  Meuser,  I.  Borovok,  M.  Betzler,  and  H.  Schrempf.  2000.  Pore 
mutations  affecting  tetrameric  assembly  and  functioning  of  the  potas¬ 
sium  channel  KcsA  from  Streptomyces  lividans.  FEBS  Lett.  472:83-87. 

Stampe,  P.,  J.  Arreola,  P.  Perez-Cornejo,  and  T.  Begenisich.  1998.  Non- 
independent  K+  movement  through  the  pore  in  IRK1  potassium  chan¬ 
nels.  J.  Gen.  Physiol.  112:475-484. 

Tanford,  C.,  and  R.  Roxby.  1972.  Interpretation  of  protein  titration  curves: 
application  to  lysozyme.  Biochemistry.  11:2192-2198. 

Terlau,  H.,  H.  Boccaccio,  B.  M.  Olivera,  and  F.  Conti.  1999.  The  block  of 
Shaker  K+  channels  by  K-conotoxin  PVIIA  is  state  dependent.  J.  Gen. 
Physiol.  114:125-140. 


van  Gunsteren,  W.  F.,  S.  R.  Billeter,  A.  A.  Eising,  P.  H.  Hiinenberger,  P. 
Kruger,  A.  E.  Mark,  W.  R.  P.  Scott,  and  I.  G.  Tironi.  1996.  Biomolecular 
Simulation:  The  GROMOS96  Manual  and  User  Guide,  vdf  Hochschul- 
verlag,  Ziirich. 

Wall,  M.  E.,  S.  Subramaniam,  and  G.  N.  Phillips,  Jr.  1999.  Protein  structure 
determination  using  a  database  of  interatomic  distance  probabilities. 
Protein  Sci.  8:2720-2727. 

Zhou,  M.,  J.  H.  Morais-Cabral,  S.  Mann,  and  R.  Mackinnon.  2001.  Potas¬ 
sium  channel  receptor  site  for  the  inactivation  gate  and  quaternary  amine 
inhibitors.  Nature.  411:657-661. 


39 


Biophysical  Journal  81(5)  2473-2483 


APPENDIX  B 

Three-Dimensional  Continuum  Simulations  of  Ion  Transport  Through  Biological 
Ion  Channels:  Effect  of  Charge  Distribution  in  the  Constriction  Region  of  Porin 


40 


nn  Journal  of  Computational  Electronics  1:  335-340,  2002 

©  2002  Kluwer  Academic  Publishers.  Manufactured  in  The  Netherlands. 


Three-Dimensional  Continuum  Simulations  of  Ion  Transport  Through 
Biological  Ion  Channels:  Effect  of  Charge  Distribution  in  the  Constriction 

Region  of  Porin 

T.A.  VAN  DER  STRAATEN,  J.  TANG  AND  R.S.  EISENBERG 
Department  of  Molecular  Biophysics  and  Physiology,  Rush  Medical  College,  1750  W.  Harrison  St., 

Chicago,  IL  60612,  USA 

trudyv@uiuc.edu 

jtang@msh.edu 

beisenbe@msh.edu 


U.  RAVAIOLI  AND  N.R.  ALURU 

Beckman  Institute  for  Advanced  Science  and  Technology,  University  of  Illinois  at  Urbana-Champaign, 

405  N.  Mathews  Ave.,  Urbana,  IL  61801,  USA 

ravaioli  @  uiuc  .edu 
alum@uiuc.edu 


Abstract.  Drift-diffusion  models  are  useful  for  studying  ion  transport  in  open  protein  channel  systems  over  time 
scales  that  cannot  be  resolved  practically  by  detailed  molecular  dynamics  or  quantum  approaches.  Water  is  treated 
as  a  uniform  background  medium  with  a  specific  dielectric  constant  and  macroscopic  current  flow  is  resolved  by 
assigning  an  appropriate  mobility  and  diffusivity  to  each  ionic  species.  The  solution  of  Poisson’s  equation  over  the 
entire  domain  provides  a  simple  way  to  include  external  boundary  conditions  and  image  force  effects  at  dielectric 
discontinuities.  Here  we  present  a  3-D  drift-diffusion  model  of  ion  (K+  and  Cl-)  permeation  through  the  porin 
channel  ompF,  and  its  mutant  G119D,  implemented  using  the  computational  platform  PROPHET. 

Keywords:  ion  channels,  Poisson-Nernst-Planck,  transport  simulation,  ompF  porin,  nanotechnology 


1.  Introduction 

Ion  channels  are  a  class  of  proteins  that  form 
nanoscopic  aqueous  tunnels  to  control  the  passage  of 
ions  through  the  otherwise  almost  impermeable  mem¬ 
branes  of  all  biological  cells.  Each  channel  consists  of  a 
chain  of  amino  acids  carrying  a  strong  and  rapidly  vary¬ 
ing  localized  charge.  From  a  biological  point  of  view, 
ion  channels  maintain  the  correct  internal  ion  composi¬ 
tion  that  is  crucial  to  cell  survival  and  function.  They  di¬ 
rectly  control  electrical  signaling  in  the  nervous  system, 
muscle  contraction,  and  the  delivery  of  many  clinical 
drugs  (Hille  1992).  There  are  many  types  of  ion  chan¬ 
nels,  each  with  a  specialized  function.  Some  channels 


have  the  ability  to  selectively  transmit  or  block  a  partic¬ 
ular  ion  species  and  many  exhibit  switching  properties 
similar  to  nanoscale  electronic  devices. 

ompF  is  a  trimeric  porin  channel  that  resides  in  the 
outer  membrane  of  the  e-coli  bacterium.  Its  well-known 
and  extremely  stable  molecular  structure  make  it  a 
good  choice  for  experimental  and  computational  stud¬ 
ies.  Each  monomer  appears  to  operate  independently, 
making  it  a  possible  candidate  for  multilevel  logic  de¬ 
vices.  ompF  carries  a  net  negative  charge  of  approxi¬ 
mately  — 30|e|,  where  e  is  the  electron  charge,  and  is 
moderately  selective  for  cations.  The  narrow  constric¬ 
tion  region  of  each  hourglass-shaped  channel  is  highly 
charged  due  to  the  presence  of  three  positively  charged 
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Figure  1.  Molecular  structure  of  ompF,  projected  along  the  length  of  the  channel,  showing  the  three-fold  symmetry  of  the  trimer.  The  five 
ionized  amino  acids  in  the  constriction  region  of  each  pore  are  highlighted. 


(R42,  R82,  R132)  and  two  negatively  charged  ( D113 , 
El  17)  amino  acids  (Phale  et  al.  2001),  as  shown  in 
Fig.  1 .  The  arrangement  of  these  charges  gives  rise  to 
a  very  strong  electric  field  parallel  to  the  plane  of  the 
membrane,  which  is  believed  to  govern  the  behavior  of 
the  open  channel.  ompF  can  be  mutated  by  replacing  or 
deleting  one  or  more  of  the  amino  acids,  thus  altering 
the  charge  distribution  along  the  channel.  Engineer¬ 
ing  channels  with  specific  conductances  and  selectiv- 
ities  is  thus  conceivable.  The  mutant  G119D ,  formed 
by  replacing  an  uncharged  amino  acid  (glycine)  with  a 
negatively  charged  aspartate,  has  a  structure  almost  un¬ 
changed  from  that  of  ompF  yet  its  conductance  is  mea¬ 
sured  to  be  between  15%  and  50%  lower  than  ompF, 
depending  on  the  concentration  of  the  salt  solution  in 
which  it  is  immmersed.  Though  the  overall  structure  of 
the  porin  is  not  compromised  by  this  mutation  the  con¬ 
striction  region  becomes  significantly  narrower  due  to 
the  bulk  and  orientation  of  the  aspartate.  The  observed 
reduction  in  channel  conductance  is  believed  to  be  a 


result  of  changes  in  both  the  steric  and  electrostatic 
environment  of  the  pore  constriction. 

Here  we  describe  a  three-dimensional  (3-D)  drift- 
diffusion  model  of  ion  (K+  and  CD)  permeation 
through  ompF  and  its  mutant  G119D,  implemented 
using  the  computational  platform  PROPHET  (http:// 
www-tcad.stanford.edu/).  While  continuum  models 
sacrifice  resolution  of  molecular  detail  they  can  be 
used  to  compute  macroscopic  current  with  a  reason¬ 
able  amount  of  computational  effort  and  have  been 
found  to  describe  permeation  through  ion  channels  sur¬ 
prisingly  well  (Hollerbach  et  al.  2000).  Using  this  ap¬ 
proach  we  investigate  the  effect  of  the  charge  distri¬ 
bution  in  the  pore  constriction  on  the  open  channel 
conductance  and  selectivity.  In  Section  2  we  describe 
our  model;  in  Section  3  we  present  the  results,  in  par¬ 
ticular  we  compare  the  current-voltage  characteristics 
and  selectivities  of  G119D  and  ompF.  Section  4  con¬ 
cludes  with  a  brief  discussion  of  this  work  and  future 
plans. 


42 


3-D  Continuum  Simulations  of  Ion  Transport  337 


2.  Poisson-Drift-Diffusion  Model 

The  channel  system  studied  here  is  comprised  of  a 
porin  trimer  in  situ  in  a  cell  membrane,  immersed  in  an 
aqueous  bath  of  KC1.  Experimentally  it  is  possible  to 
maintain  different  salt  concentrations  on  either  side  of 
the  membrane  (hereafter  referred  to  as  and  Cr, 
Electrodes  are  immersed  in  the  baths  to  apply  a  fixed 
bias  across  the  channel/membrane  system.  The  elec¬ 
trostatic  environment  of  the  channel  is  determined  by 

(i)  mobile  ions  permeating  through  the  open  channel, 

(ii)  fixed  (permanent)  charges  that  reside  on  the  pro¬ 
tein  itself,  and  (iii)  the  charges  in  the  aqueous  baths 
and  on  the  electrodes.  The  local  electrostatic  potential 
c p  is  related  to  all  the  charges  in  the  system  by  Poisson’s 
equation, 

V  •  (eV<?)  =  -{Pfaed  +  P+  +  P-)  (1) 

where  e  is  the  dielectric  constant,  and  Pfixed,  p+  and 
p  are  the  charge  densities  per  unit  volume  of  fixed 
charges  residing  on  the  protein,  K+  ions  and  Cl-  ions, 
respectively.  Current  flow  j±  is  described  by  the  drift- 
diffusion  equation 

j±  =  -(p±p±V<p  ±  D± V p± )  (2) 

where  n±  and  D±  are,  respectively,  the  mobilities 
and  diffusivities  of  the  ionic  species.  Conservation  of 
charge  dictates 

-*  dp± 

v  •  j±  +  -p  =  S±  (3) 

at 

where  S±  can  be  any  form  describing  the  details  of  ion 
binding  and  other  chemical  phenomena  that  populate 
or  deplete  the  ion  densities.  In  the  present  model  we  do 
not  consider  such  phenomena  and  set  S±  —  0.  We  seek 
a  steady-state  solution  for  tp,  p+  and  p_that  simultane¬ 
ously  satisfies  Eqs.  ( 1)— (3)  subject  to  specific  Dirichlet 
boundary  conditions  at  the  electrodes  (C/f/,,  C  right  and 

F bias)  • 

The  model  described  above  was  implemented  us¬ 
ing  the  PROPHET  simulator,  a  computational  plat¬ 
form  originally  developed  at  Lucent  Technologies 
and  currently  being  extended  at  Stanford  University. 
The  PROPHET  simulator  uses  the  “dial-an-operator” 
methodology  to  construct  systems  of  partial  differen¬ 
tial  equations  by  combining  existing  differential  opera¬ 
tors  described  using  a  scripting  syntax.  Equations  (1)- 
(3)  and  the  boundary  conditions  are  constructed  using 


existing  PROPHET  operators.  After  the  channel  ge¬ 
ometry  is  defined  on  a  customized  mesh  provided  by 
the  user,  the  equations  are  discretized  and  the  lin¬ 
ear  system  is  solved  using  iterative  methods.  Phys¬ 
ical  parameters  are  assigned  to  the  different  regions 
of  the  domain  at  run-time.  For  the  results  presented 
here  we  use  s  =  80,  20  and  2,  for  the  aqueous,  protein 
and  membrane  regions,  respectively.  Diffusivities  were 
chosen  to  fit  measured  current-voltage  (/-VO  curves, 
assuming  D+  =  />_,  and  mobilities  were  assigned  us¬ 
ing  the  Einstein  relations,  reducing  the  number  of  ad¬ 
justable  parameters  to  one. 

The  distribution  of  charge  on  each  amino  acid  is 
modelled  by  associating  a  fractional  point  charge  (in 
units  of  \e\)  to  each  atom.  Values  for  these  charges  in 
neutral  solutions  are  calculated  using  ab  initio  quantum 
chemistry  codes.  For  this  work  we  have  used  the  tabu¬ 
lated  OPLS  partial  charges,  a  force  field  used  widely  in 
molecular  dynamics  simulations  of  protein  molecules 
(Jorgensen  and  Tirado-Rives  1988). 

3.  Results 

As  mentioned  in  the  previous  section,  the  distribution  of 
permanent  charge  around  each  amino  acid  is  obtained 
from  quantum  chemistry  calculations  that  treat  each 
amino  acid  as  a  free  entity  in  solution.  The  dielectric  en¬ 
vironment  and  proximity  of  neighboring  ionized  amino 
acids  in  the  folded  protein  is  expected  to  alter  the  charge 
distribution.  Electrostatic  calculations  by  Karshikoff 
et  al.  (1994)  suggest  that  at  neutral  pH  the  charges  on 
amino  acids  R82,  located  inside  the  constriction  region, 
and  R 1 67 ,  located  near  the  mouth  of  the  channel,  should 
be  scaled  down  from  +1  to  zero  (Dutzler  private  com¬ 
munication).  More  recent  calculations  (Schirmer  and 
Phale  1999,  Varma  and  Jakobsson  private  communica¬ 
tion)  however,  indicate  that  both  R82  and  R167  should 
in  fact  be  fully  ionized  at  neutral  pH.  Since  both  R82 
and  R 1 67  line  the  pore  of  the  channel  it  is  important  that 
these  charges  be  properly  represented.  Figure  2  shows 
I-V  curves  computed  for  symmetric  bath  concentra¬ 
tions  of  (a)  100  mM  KC1  and  (b)  1  M  KC1,  with  both 
representations  of  the  charge  distribution:  R82,  R167 
uncharged — solid  symbols;  R82,  R167  fully  ionized — 
open  symbols.  Figure  2(c)  shows  the  computed  I-V 
curve  for  a  system  with  100  mM  KC1  on  one  side  of 
the  membrane  and  1  M  KC1  on  the  other,  assuming 
R82  and  R167  are  uncharged.  A  diffusivity  in  the  range 
(0.8-1.02)  x  10~5  cm2  s“  1  provides  a  satisfactory  fit  to 
the  experimental  data.  Interestingly,  the  extra  positive 
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(a) 


(b) 


vbias<mV> 

(c) 


Figure  2.  Comparison  of  drift-diffusion  model  with  experimental 
I-V curves  for  ompF in  (a)  100  mM  (b)  1  M  and  (c)  100  mM  :  1  M  KC1. 
Simulations  run  with  both  R82  and  R167  fully  ionized  are  indicated 
with  open  circles. 

charge  on  the  protein  does  not  make  a  significant  dif¬ 
ference  to  the  open  channel  conductance.  However  the 
selectivity  of  the  ompF  for  cations,  defined  here  as  the 
fraction  of  total  current  carried  by  cations,  is  reduced 
from  66%  to  62%  at  1  M  KC1,  and  88%  to  84%  at 


Figure  3.  Comparison  of  drift-diffusion  model  with  measured 
I-V  curves  for  G119D  in  100  mM  and  1  M  KC1.  Measured  data  are 
indicated  with  solid  (100  mM)  and  dashed  (1  M)  lines,  the  simulation 
data  are  indicated  with  symbols. 

100  mM  KC1,  which  consistent  with  the  addition  of 
+6\e  \  making  the  channel  less  electrostatically  accessi¬ 
ble  to  cations.  These  values  are  commensurate  with  the 
value  77%  measured  for  100  mM  -  1  M  KC1  (Schirmer 
and  Phale  1999). 

G119D  is  created  by  replacing  the  neutral  glycine, 
G119,  tucked  between  the  aspartate  Dll 3  and  gluta¬ 
mate  El  17  (see  Fig.  1),  both  negatively  charged,  with 
an  aspartate  D119,  also  negatively  charged.  In  addi¬ 
tion  to  increasing  the  negative  charge  the  bulkier  aspar¬ 
tate  also  reduces  the  volume  available  for  ions  to  pass 
through  the  constriction  zone.  Figure  3  shows  the  mea¬ 
sured  and  computed  /-^characteristics  for  G119D.  At 
100  mM  KC1  the  conductance  is  about  15%  lower  than 
that  of  ompF,  while  at  1  M  the  conductance  of  G1 1 9D  is 
about  40%  lower.  A  diffusivity  in  the  range  (6.0-6. 8)  x 
10~6  cm2  s^1  is  able  to  fit  a  fairly  large  range  of  experi¬ 
mental  data.  The  simulations  predict  G1 1 9D  to  be  more 
cation  selective  than  ompF  (69%  at  1  M  and  90%  at  100 
mM  KC1),  as  would  be  expected  from  the  addition  of 
an  extra  charge  in  the  constriction  region. 

The  mutation  introduces  both  steric  and  electrostatic 
changes  in  the  constriction  region  of  the  channel.  How¬ 
ever,  the  relative  contribution  of  each  change  to  the  re¬ 
duced  conductance  observed  in  G119D  is  not  known. 
Understanding  the  relationship  between  structure  and 
function  of  ion  channels  is  crucial  if  new  channels  are 
to  be  designed  with  specific  properties  and  behaviors. 
Simulations  can  perhaps  find  their  greatest  use  in  pro¬ 
viding  insight  that  cannot  be  measured  or  inferred  from 
experiment.  As  an  example  of  this  we  have  attempted  to 
separate  the  steric  and  electrostatic  effects  in  G119D  by 
scaling  down  the  charge  distribution  on  the  aspartate 
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to  zero.  This  creates  a  channel  (albeit,  a  virtual  one) 
that  should  be  electrostatically  similar  to  ompF  but 
structurally  the  same  as  G119D.  Simulations  at  1  M 
KC1  reveal  that  the  conductance  (shown  in  Fig.  3)  of 
such  a  channel,  were  it  to  exist,  is  almost  the  same  as 
G119D  despite  the  difference  in  the  charge  constella¬ 
tion  in  the  constriction  region.  The  selectivity,  on  the 
other  hand,  reduces  to  56%.  This  raises  the  question 
of  whether  the  lowered  conductance  of  G1 1 9D  is  more 
the  result  of  a  narrowing  of  the  pore  constriction  than 
any  electrostatic  changes.  One  issue  that  needs  further 
probing  is  how  the  orientation  of  the  key  amino  acids  in 
constriction  region  depends  on  the  local  field.  If  it  were 
possible  to  remove  the  charge  on  the  aspartate  would 
the  charged  amino  acids  reorient  in  such  a  way  that 
the  channel  cross-section  is  appreciably  altered?  De¬ 
tailed  molecular  dynamics  simulations  are  required  to 
answer  these  types  of  questions.  An  interesting  muta¬ 
tion  to  consider  experimentally  would  be  G1 1 9 /V  where 
the  uncharged  glycine  is  replaced  with  an  asparagine,  to 
emulate  the  bulk  of  the  aspartate  without  the  charge.  If 
the  above  conjecture  is  correct  then  the  conductance  of 
G1 1 9N  should  be  similar  to  G1 1 9D  and  its  selectivity 
similar  to  ompF. 

Ion  occupancies,  defined  as  the  integral  of  the  ion 
density  over  a  specified  volume,  have  been  calculated 
for  ompF  and  G1 1 9D,  for  a  range  of  salt  concentrations 
and  applied  bias.  The  latter  has  very  little  effect  on 
the  predicted  ion  occupancies,  since  the  effective  fixed 
charge  lining  the  pore  is  so  strong.  Figure  4  compares 
the  cation  and  anion  occupancies  for  the  four  channel 
scenarios  discussed  above,  computed  for  the  conditions 
of  zero  applied  potential  and  1  M  KC1.  In  all  four  cases, 
occupation  of  the  constriction  region  is  more  favorable 
for  cations,  and  the  predicted  changes  in  occupancy 


Figure  4.  Effect  of  adding  or  deleting  charges  in  the  channel 
constriction  on  the  cation  and  anion  occupancies. 


upon  mutation  are  consistent  with  the  changes  in  the 
net  charge  in  the  pore  constriction. 

Conclusion 

The  success  of  the  drift-diffusion  theory  in  predicting 
measured  TV  curves  of  ompF  and  G119D  is  reasonably 
good  considering  that  we  have  used  only  one  adjustable 
parameter  (diffusivity)  to  simulate  ion  permeation  in  a 
wide  range  of  experimental  conditions.  The  calibrated 
diffusivities  agree  to  within  a  factor  of  3  with  the  exper¬ 
imentally  measured  values  in  bulk  salt  solution.  Theo¬ 
retical  values  for  the  ion  diffusivities  inside  the  porin 
channel,  as  calculated  from  molecular  dynamics  simu¬ 
lations,  also  agree  with  the  calibrated  values  to  within 
20%.  Clearly  a  uniform  diffusivity  cannot  capture  all 
the  physical  effects  affecting  ion  transport  in  an  inher¬ 
ently  three-dimensional  restricted  and  highly  charged 
volume.  The  development  of  a  diffusivity  model  for 
each  ionic  species  that  includes  a  dependence  on  posi¬ 
tion  (or  local  field)  and/or  ion  density  is  the  subject  of 
ongoing  work. 
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